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A parameter vector (N components a 2 , . . . , a N ) 

ith parameter 
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MN(i) dependency index associated with the parameter a^ (FORTRAN 
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n-L index of the ith vector of the optimal basis and of the 

corresponding parameter 
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P number of outputs 
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s D separation threshold 

time instants 

U input vector (forcing function) 
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X state vector 

Y measured output vector (M components y(ti), y(t 2 ) , . . 

y(t M )) 

Yq projection of Y upon % 


y(or y(t)) measured time history of an output 
e angular threshold 

2 ft model parameter subspace (hypersurface in output space Em) 

( ) time derivative of ( ) 

( ) estimated value of ( ) 

( ) transposed of the matrix ( ) 

( ) 1 inverse of the matrix ( ) 


b 

C 

E 

g 
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h W 


Symbols Used in the Illustrative Example 
bearing offset * mass M of the platform (three components 



control matrix (diagonal 3 x 3 of elements C ls C 2 , C 3 ) 

voltage applied to the DC motor 

gravity vector (3 components) 

body angular momentum vector (3 components) 

total angular momentum of the reaction wheels (3 components 
h Wl’ h W 2 ’ h W3^ 
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current in the DC motor 

inertia tensor (inertias , Jx, Jy, Jz; products, Jxy, Jxz, Jyz) 
inertia tensor of the wheel-rotor combination (diagonal) 
back electromotive force constant 
damping constant 

torque constant of the DC motor 
DC gains of the filters 
resistance of the wiring 
skew matrix defining a cross product 
total torque applied to the body 
gravity torque 

total torque due to the reaction wheels 

feedback matrix (diagonal 3 x 3 of elements 04 , a^, a 3 ) 
time constant of the first order filters 
body angular velocity vector (3 components p, q, and r) 
angular velocity of the wheel 

Symbols Used in the Appendixes 
jth component of the vector D^ ER 
product of the i first pivots 
unit vector parallel to 8Y/3a ^ 

differential of a^ or error bound on the parameter a^ 

expectation value of ( ) 

unit vector orthogonal to some space 

linearly close with respect to the threshold s Q 

value of the i th pivotal element 

separation of some vector D from some set B 

vii 



■mi ii iiiii mu ■ ■ ■■ ■ ■■■■■! ■■ ■ 

s i 

separation of some vector from some set 


V 

some vector of Ek-1 


V (B) 

volume of the hyperparallelepiped spanned by the vectors of a 
set B 

Z 

function of z ^ 


^1 j z 2 j ■ 

coordinates in the space ^ ' 

K-l 


SY 

error in Y 


aY o 

error in Y Q 


e K 

part of 3Y/3a K orthogonal to Ej^ 


CM O 
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square of the error bound in Y 0 


*i 

component of the vector 3Y/3a on the vector 

aY/aai 

y i 

component of the vector on the vector 
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l 

component of the vector D K on the vector 
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in another 

U 

set union 


c 

set inclusion 
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set membership 


V 

universal quantifier (for any) 
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NONLINEAR SYSTEMS IDENTIFICATION IN PRESENCE OF NONUNIQUENESS 

Jean-Noel Aubrun* 


Ames Research Center 


SUMMARY 


This report considers the problem of identifying a set of parameters that 
will match the input-output of a mathematical model to that of a physical sys- 
tem. Except for particular cases, there has been no practical method to 
determine if such parameter values are unique. A general process of parameter 
identification is described that addresses this problem and uses geometrical 
concepts in terms of which the nonuniqueness problem can easily be defined. 

A digital computer algorithm is developed that analyzes the structure of 
a space defined by the parameter sensitivity functions and the output data set 
The algorithm deduces an optimal set of parameters to be uniquely identified, 
determines the relationships between dependent parameters, and specifies which 
parameters can be obtained with a priori knowledge of others. It does not 
require canonical or linear equations for the model but maintains the physical 
identity of the parameters. A corresponding FORTRAN IV program has been 
written. 

This technique is illustrated by the identification of the Ames three- 
degrees-of-freedom satellite simulator. Examples of nonuniqueness were found 
and analyzed successfully by the algorithm, demonstrating its ability to cope 
with strongly nonlinear cases. 


INTRODUCTION 


The relationship between the input and the output of a physical system or 
plant may be described in most situations by a mathematical relationship or 
model. This paper is concerned with models that depend on a finite set of 
parameters and with the identification of these parameters from measurements 
of the plant input and output. 

General theories of parameter identification have been established (e.g., 
refs. 1 and 2), and specific applications have been made (e.g., refs. 3-7) 
primarily to linear systems for which a thorough mathematical analysis can be 
performed. For this discussion, a distinction should be made between two 
types of parameter identification. In the first type, the goal is to obtain a 
mathematical representation of the plant from which it is possible to dupli- 
cate input-output sets of measurements regardless of the actual physics or 
dynamics involved in the plant. This is the concern of realization theory. 

The algorithm of Ho (ref. 7) and later developments introduced by 

^National Research Council Postdoctoral Associate in residence at Ames 
Research Center. 



Kalman (ref. 4) are typical of this approach. In the second type the form of 
the model is obtained first from an analysis of the physical processes occur- 
ring in the plant. The parameters in this case have a direct physical meaning 
and the goal is to obtain their actual values from input-output measurements. 

If these values are known exactly, plant and model will exhibit the same input- 
output properties. However, the converse is not true although often implicitly 
assumed (i.e., the parameter values that match model and plant input-output 
may not be unique) . 

This report is mainly concerned with the second type of parameter 
identification where it is essential to recognize the occurrence of nonunique 
solutions for the parameter values. This nonuniqueness may come from param- 
eter redundancy or insufficient information in the measurements or both. 

While it is possible in linear systems to transform the original parameters 
into a uniquely identifiable reduced set (canonical parameters) (ref. 8), 
there is no general canonical solution for nonlinear systems. Concerning the 
second cause of nonuniqueness, theoretical criteria are available from estima- 
tion theory for linear systems, but because they involve idealized arithmetic 
operations (such as identically vanishing determinants), they may fail com- 
pletely when numerical computation or noisy measurements are used. 

To solve the nonuniqueness problem, this report first presents a general 
procedure of parameter identification in a more direct and practical way than 
has been done previously. This procedure is used as a framework for analyzing 
the mechanics of parameter identification on a practical level and allows a 
formulation of the nonuniqueness problem that takes into account the uncer- 
tainties in real measurements and computations. These uncertainties are 
expressed in the relaxation of strict mathematical concepts to the benefit of 
approximate ones. It is realized, for instance, that there is no such thing 
as a singular matrix with a digital computer and the classical "linear depen- 
dence" is replaced by a concept of "linear closeness." 

A digital computer algorithm is developed that: 

1. Detects the existence of nonuniqueness. 

2. Determines which sets of parameters are not uniquely defined and 
which are independent. 

3. Determines which parameters should have values specified to obtain a 
correct answer for the others. 

4. Optionally, resets the main identification program with a reduced 
number of parameters that can be uniquely identified, the others being left at 
a constant value. 

This computation technique is applied to data recorded from the Ames 
three-degrees-of-freedom satellite simulator. Identification of the inertia 
tensor and some linear and nonlinear control elements is performed to provide 
a typical illustration of the nonuniqueness problem. 
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BASIC PARAMETER IDENTIFICATION PROCEDURE 


The Plant 

The direct measurements of a physical system (plant) may be considered as 
a set of quantities that are recorded at different instants of time. It is 
convenient here to partition the set into input and output sets. For dynamic 
systems, the input is usually defined as some forcing functions that modify 
the state of the system. This restriction is unnecessary for the present 
approach to the identification problem. Here it will be required only that 
the plant be some kind of operator connecting input and output. In the case 
of an airplane, for instance, longitudinal stick position could be the input 
and pitch rate the output. However, depending upon the type of identification 
to be performed, one could also use longitudinal stick position and pitch rate 
as input, and pitch acceleration as output. 


The Model 

To define the plant completely one must know all the output sets that 
correspond to all the possible input sets. However, this is not very practi- 
cal. Instead one tries to determine what is called a "model," that is some 
process by which, for any given set of input values, one is able to obtain the 
corresponding set of output values. It will never be possible to do so 
exactly, first because of the actual complexity of the plant and because a 
finite number of digits have to be used in computation. Hence the model is 
only an approximation that describes the plant to the extent that the sets of 
input and output values of the model are close in some sense (such as least 
squares) to those of the plant. The approximation, however, will focus on the 
important properties of the plant, so that its behavior can be predicted or 
modified in some desired way. 

The model in itself depends upon quantities called the parameters. In 
this discussion we will consider only models whose parameters are invariant 
with respect to time and of finite number. For instance the mathematical 
model of an aircraft may be a system of differential equations describing its 
dynamics. The coefficients in these equations, such as stability derivatives, 
inertias, and control gains along with the initial conditions, are normally 
considered the parameters. 


The Identification 

The objective of the identification process is to determine, for a given 
model form, the values of the parameters that will make the model behave like 
the plant. A general process of identification is sketched on figure 1. If 
the same input is applied to both model and plant, the best values for the 
parameters are those that minimize (with respect to some criterion) the 
difference between model and plant outputs. Therefore the identification 
algorithm has to compare these outputs and adjust the parameter values until 
this minimum is reached. 
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It is usual, in parameter identification terminology, to distinguish 
between "equation error" and "output error" methods (ref. 6). It should be 
pointed out that both of them can be cast into the scheme of figure 1, 
provided a judicious choice of input and output is made. Let us consider for 
instance a plant defined by the system of differential equations: 

X = F (X, U) (1) 

where X is the state vector and U the control. An output error method 
will consider U as the input to a model defined by 

X = F (X, U) 

and compare the output X of the model to the measured output X of the 
plant to obtain the parameters in F. In an equation error method, U and X 
can be considered as inputs to a model defined by 

X = F (X, U) 

while X is now the model output to be compared to the plant output X. 

In the first case the model is differential, in the second it is 
algebraic, but in both the only information concerning the parameters is to be 
found in the comparison between measured and computed outputs. 

Once the form of the mathematical model is specified, two important 
questions must be asked: 

1. If it is possible to match some particular set of input-output values 
of the system to that of the model, using certain values for the parameters, 
and is it possible to match every other set with the same parameter values? 

2. Is it possible to match the same input-output set with more than one 
set of values for the parameters? 

The first question relates to the problem of modeling; the second expresses 
the problem of uniqueness. A general way to perform the identification must 
be studied to understand where the uniqueness problem comes from. 

Consider first the case of a single output, say y. Measurements of the 
time history y(t) have been made, corresponding to a known input. In this 
analysis, digital computation will be considered; hence the output set con- 
sists of M discrete values of y corresponding to the instants ti, t 2 , 

. . ., tjq, as represented on the top curve of figure 2(a), dotted at the 
measurement points. The model depends upon the N parameters of unknown 
values, ai, a 2 , . . ., To begin with, these parameters are given some 

assumed values, purely arbitrary or based upon a best guess. Since the input 
is known, a computed time history y(t) may be obtained from the model, 
particularly the M discrete values corresponding to the same instants ti, 
t 2 , . . ., tM* In general, the two time histories of y and y are different 
as shown also on figure 2(a), and the problem is to find how the parameters 
must be adjusted to make these time histories as similar as possible. 
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An intuitive idea is to change the value of one parameter, say aj, by a 
small amount dai , and observe what happens to the time history of y. A 
slightly different time history y + dy x is obtained (fig. 2(a)), and even- 
tually some idea of the effect of the parameter and how much it should be 

changed could be gained. The same can be done with a 2 , obtaining another 
time history with different characteristics, and so on with the N parameters 
of the model. By this trial and error procedure, one might be successful in 
matching the two time histories. This technique, called analog matching, has 
been fairly successful in the past. It is obvious, however, that it becomes 
impractical if the parameters are numerous, although the human mind takes 
definite advantage of its pattern recognition ability. 

This process may be rationalized by considering the M values of y as 
the components of a vector Y,^and the M corresponding values of y, as the 
components of another vector, Y (fig. 2(b)). Both Y and Y belong to an 
M-dimensional space and there is a one-to-one correspondence between a time 
history and its representative vector. To say that the time histories are 
different is to say that the vector Y and Y are different. Their difference 
may be expressed directly by the geometric difference Y - Y, and we define 
therefore an M component error vector as 

ER = Y - Y 

It follows that the discrete time histories will match iff ER is zero. 
Consider now what happens when the value of is changed. The change in 

the time history noted previously has now a very precise meaning, which is the 
change d?i in the vector Y that becomes equal to 

AAA 

Y x = Y + dY ! 

A 

This M component vector dY* describes unambiguously how the time history 
has been modified by the increment da^ of the parameter a^ . If ds-i is 
small enough, one may write 



and all the local information concerning the effect of a^ upon the time 
history is contained in the vector 3Y/3ai- Therefore, the first step in the 
identification will be to perturb each parameter a^ individually and compute 
the corresponding derivatives (M components vectors) 3t/3a^ by 

9 ^ = (Yj - 
3a^ da^ 

It has just been demonstrated that if the values of the parameters are 
changed, it is possible to change the vector Y in this M-dimensional space. 
It would be desirable to change it in such a way that it will coincide with Y 
(fig. 3). If simultaneous ^increments da^ , da 2 , . . da^j are given to the 
corresponding parameters, Y will change by an amount 
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d? -Et(|y dM 


i = 1, N 


( 2 ) 


This relation may be expressed conveniently in matrix form by considering each 
parameter as one component of an N-dimensional vector A, and each vector 
9f/8a^ as a column of an M x N matrix DER. Then equation (2) is 
equivalent to 

dY = DER dA (3) 

with 


DER 


(ti)/9ai 

9y(t 2 )/3ai 


>3y(tM)/ 3a l 


3y(tl)/3 a 2 

9y(t 2 )/9a 2 


9yCtM)/3 a 2 


9y(ti)/9a N ' 

9y(t 2 )/9a N 


9y(tM)/3 a N, 


Finally a dA is needed such that the corresponding change in Y just 
cancels the error ER, that is, such that 


ER + DER dA = 0 (4) 

The left side represents the error after the parameter vector has been 
changed by dA, resulting in the change dY in Y as shown on figure 3. 
Usually there are more measurements than parameters (M >> N) so that equa- 
tion (4) cannot be solved directly. However, a least squares solution may be 
obtained by minimizing the quantity (ER + DER dA) 2 with respect to dA. A 
well-known expression for the solution is given by 

dA = -(DER t DER) -1 DER T ER (5) 


This solution is valid, of course, provided the inverse of the square N x N 
matrix (DER^ DER) exists. 

This type of equation is found in many identification techniques (refs. 1, 
2, 3, 6). It represents a pseudoinversion of ^DER and has an interesting 
geometrical meaning. As has been said, Y and Y belong to an M-dimensional 
space (fig. 4) and it is possible to "move" ? in this space by changing the 
value of the parameters. But, since there are only N parameters (N < M) , 
they define at most N independent directions; that is, Y is, in fact, con- 
strained to stay in an N-dimensional subspace of the M-dimensional 

space. If there were no noise in the measurements, if the model were perfect, 
and if it were possible to compute everything exactly, the vector Y would 
also belong to this subspace 2^ and Y and Y could be made to coincide 
exactly, cancelling completely the error ER and therefore obtaining a per- 
fect match of the time histories. This situation corresponds algebraically to 
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compatible equations in the system expressed by equation (4). Unfortunately 
the conditions expressed above are never met in practice, so Y is always 
outside this subspace 2^ and cannot be reached by Y. Therefore the best 
tack is to reach the orthogonal projection Y 0 of Y upon Sjsj. The distance 
[ (Y - Y ) 2 ] 1 / 2 will then be minimal so that this can be interpreted as a least 
squares solution. Equation (5) expresses nothing but that the real M target M 
of Y is not Y but Y 0 . Indeed from equations (3) and (5) 

dY = DER dA = -DER(DER T DER)” 1 DER T ER = -ER 0 

It can be shown that ER 0 is just the orthogonal projection of ER on % 
by verifying that (ER - ER 0 ) is orthogonal to all the vectors 3Y/3Aj_ which 
constitute a local basis in 2^; that is, since these vectors are contained in 
the matrix DER, 

DER t (ER - ER 0 ) = DER T ER - DER T DER(DER T DER)” 1 DER T ER = 0 

It must-be emphasized, however, that equation (5) is not the only way to 
operate on Y. There are many different techniques for obtaining a least 
squares solution of equation (4) and they may not all require an actual matrix 
inversion as in equation (5). But it is very important to note that the^pro- 
jection property is independent of any technique; 1 whatever the method, Y 
will always be constrained in some ^N-dimensional subspace 2^ and the param- 
eters adjusted in such a way that Y can reach some target Y 0 in 2^. Also 
independent of the identification technique itself is the fact that any change 
dY in the output of the model has to be related to the change dA in the 
parameter vector by 

dY = DER dA 

This expression defines the local properties of the space 2^. The extremity 
of Y is therefore constrained on an N-dimensional hypersurface 2^. The 
local hyperplane E^ tangent to 2^ is spanned by the column vectors of DER 
(i.e., the N vectors 9?/8Ai). When DER is constant with respect to the 
parameters, 2^ is also a hyperplane and coincides with Ejq everywhere. When 
it is not, 2 n is a curved hypersurface, and equation (5) will not move ¥ 
toward Y 0 but toward the projection Y 0 of Y on Ej\j. Because Y 0 is 
only ^an^approximation for Y Q> equation (5) must be iterated. The evolution 
of Y, Y 0 , and En during this iterative process is represented schematically^ 
in figure 5; Y has to move^along the curve 2^ and the tangent to 2^ at Y 
is En. At the beginning Y is Y(0). Equation (5) approximates Y 0 by the 
projection Y 0 (l) of Y upon En(0) x A^new value tf(l) is obtained for Y 
at the first iteration and the arc (Y(O)Y(l)) on 2^ is approximately equal 
to the distance Y(0)Y o (l) on En( 0) if ?(0) were close to Y 0 . Starting 
now with Y(l), a better approximation of Y 0 is obtained, Y 0 (2). Finally, 
when Y(i) is close enough to Y 0 , then Ejsj(i) and 2jsj are equivalent to 


1 When weights are used in the minimization of the error (as in maximum 

likelihood estimation or in general weighted least squares methods) , these 
results remain valid provided the output vector is redefined as, say Y 1 , such 
that Y ! = WY where W is a weighting matrix. 
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compute the projection of Y, and, in general, any linear property of 2^ 
will also be found in Eft. 

Once Y 0 is reached, that is, when the time histories are matched as 
well as possible (i.e. , when ER 2 is minimum, ER 2 being the cost of the 
identification since it is a measure of how good the match is), one must ask 
if it is possible to change the parameters in such a way that the output Y 
does not change. If the answer is yes, then it is possible to find another 
set of values for the parameters that will match the time histories as well; 
hence the parameters cannot be uniquely identified. Since 2^ contains all 
the properties of the model, it is clear that this nonuniqueness has to be 
related to the properties and the structure of 2^, or to those of E^ if 
only linear properties are involved. 

Before this point is developed, a last word should be said for the 
multioutput case. If the system has P outputs, y x , y 2 , . . . , yp, each com- 
posed of M 0 measurements, then the vector Y is constructed with these 
M 0 x p quantities. The order in which these measurements are taken to form 
the components of Y is of no importance provided the same order is used for 
the model output Y. An example of the vectors and matrices involved is given 
in expression (6) where Y and Y are M-dimensional with M = M 0 x P and 
(Y - Y) will be zero iff the P time histories of the model are equal to 
the corresponding time histories of the system. With this definition of Y 
and Y, all the previous discussions and results are valid. 


Plant 

output 

Y(M) 

Model 

output 

Y(M) 

Matrix of the 

derivatives of 
DER (M x N) 

the model output 

yj(ti) 


yi (ti) 


3y i (ti)/3a 1 

3yi (ti)/3a 2 

• * 9yi(ti)/9a N 

yi(*2) 


yi ^ 2 ) 


3yi (t 2 )/3ai 

9y i C 1 2 3 / 9a 2 

. . dy i (t 2 ) / 3a^ 

yi (tM 0 ) 


yi (tM 0 ) 


9 y 1 (tM Q )/9ai 

9yi(tM 0 )/9a 2 • 

* • 9 yi( t M 0 )/ 9a N 

y2(ti) 


y2(ti) 


9y2 Cti3/3ai 

9y2C t l)/9a 2 

. . 3y2(ti)/9a N 

y2(*2) 


y2(t2) 


9^2 C 1 2~) / 3ai 

9y2(’ t 2)/ 9a 2 

* • 3y 2 (t 2 )/3a^ 

y2( t M 0 ) 


y 2 (tM 0 ) 


9y2( t M 0 )/9ai 

9y2( t M 0 )/ 9a 2 • 

• • 9y2( t M 0 )/ 9a N 

yp(ti) 


yp(ti) 


9yp (t]3/3ai 

9 yp(ti)/3a 2 

* • ^ypCt^/Ba^ 

yp(t 2 ) 


yp(*2) 


9yp ( 2 . 1 / 9a l 

9 yp(^2)/9a 2 

• • 9 yp(i2)/ 3a N 

>>(%) 


y P (tM 0 ) 


9 yp( t M 0 )/ 9a i 

9yp^ t M 0 )/ 9a 2 • 

• * 9yp(tM 0 )/ 9a N 
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THE NONUNIQUENESS PROBLEM 


To avoid any confusion with other possible theoretical definitions, we 
have to make clear that the nonuniqueness problem considered here is a practi- 
cal one. The only available information is a set of input-output values of 
the plant from which the parameters of the model have to be identified. This 
could be obtained, for instance, during a flight test and the corresponding 
piece of data is of finite length. It might happen that the identification 
using this particular data has nonuniqueness problems, while another piece of 
data taken with the same airplane leads to unique values. Therefore the non- 
uniqueness is relative to a given set of data, although it might also be an 
intrinsic property of the model, in which case it will always be found. 

When discussing the identification process it was said that the N 
parameters of the model defined locally an N-dimensional subspace E^, the 
basis vectors of this subspace being the M-component vectors 3?/9ai. In 
fact this is true only if the basis vectors are linearly independent. Other- 
wise E]sj is not truly N-dimensional and the fact that one vector of the 
basis can be expressed as a linear combination of some others means that the 
effect of one parameter change could be obtained also by a simultaneous change 
in some other parameters. Obviously there will not be a unique solution in 
this case. If equation (5) is used, the matrix (DERT DER) is singular and the 
computation might stop at this time. Unfortunately, the computer will very 
likely invert it anyway, because roundoff errors make it generally impossible 
to obtain a hard zero for the determinant. The solution obtained in this case 
might very well match the time histories because of the adaptative nature of 
the algorithm, as we have actually observed on simulated data. If other 
techniques are used that do not involve an actual inversion (as in the steep- 
est descent method), then a solution is reached anyway. Here again the non- 
uniqueness will be unnoticed and wrong conclusions will be drawn if one 
observes only the fit in the outputs. One should be careful not to judge 
things only for their outer properties, but to understand also the inner 
aspect (i.e., in this case, the real structure of the space En) . 

Consider a two-dimensional example to begin with, corresponding to a 
model with two parameters a x and a 2 . Here, En is E 2 (i.e x , a plane) and 
is shown on figure 6(a) with the basis vectors 9Y/9ax and 9Y/8a 2 . The 
"target" is Y 0 , the projection on E 2 of some M-dimensional vector Y. The 

amounts dax and da 2 by^which the parameters have to be changed are equal to 

the components of Y 0 - Y on the basis, because by definition of the 
components 

Vo - * - (It) * (4) 

and this is also equal to the change in Y according to equation (2) . 

Consider now the case of figure 6(b) where the basis vectors are parallel. 

Here Y 0 may be reached by changing either ax or a 2 , or both. The vectorial 

equation (7) collapses in one algebraic equation that of course has an 
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infinity of solutions in da^ and da 2 . The model does not really depend 
separately on ai and a 2 but on some function of them (see appendix A) . How- 
ever, if the true value of, say, a^ is known (from other sets of measurements 
for instance), it is possible to continue the identification process and 
obtain the true value of a 2 , a^ being removed from the perturbation process 
and left at its true value. These basic ideas will be extended to the 
N-dimensional case, but we have been dealing so far with a strict mathematical 
definition of linear dependence and independence. Here there is a definite 
boundary between the cases of figures 6(a) and 6(b). 

In the real world this yes-or-no situation does not exist and 
intermediate states must be considered. Because of noise in measurements, the 
actual output of the plant is not Y 0 but may be somewhere in a domain sur- 
rounding Y 0 . Assuming some kind of boundaries for this uncertainty domain, 
as represented on figure 7 by the shadowed area, there are corresponding 
uncertainties in the vectors dYj and dY 2 (shadowed parts of the axis in 
figure 7) that in turn correspond to uncertainties in da^ and da 2 . If the 
angle between the two vectors 3f/3ai and 3Y/3a 2 is decreased, these uncer- 
tainties increase. If they become, say, one order of magnitude larger than 
the parameters themselves, for all practical purpose these parameters are not 
uniquely defined. However, the values obtained for these parameters are not 
independent of each other, since they must satisfy the constraint that 
dYj + dY 2 falls inside the uncertainty domain of Y 0 . Therefore, if the 
angle between the two basis vectors is less than some threshold e, nonunique- 
ness will be experienced almost as if this angle were zero, and these vectors 
are called linearly close by analogy with the situation encountered with 
"linearly dependent" vectors. 

These ideas can be extended in a K-dimensional case and are developed in 
appendix B. Consider a space (fig. 8) spanned by the independent vec- 

tors 3?/3ai, 3Y/3a 2 , . . ., 3Y/3aj(-i, and a K th vector, 3Y/3a^. The K 
vectors 3Y/3ai , ^3Y/3a 2 , . . ., SY/Sa^.j , 3Y/3a^ are linearly close if the 
angle between 3Y/8a]< and its projection on the space E]<_i is smaller than 
some threshold e. It is also convenient to consider the distance between the 

A A 

extremity of 3Y/3aj( and the space Ej^j. If 3Y/3a^ is normalized to unity, 
this distance, s, is a measure of the "separation" between 3Y/3a« and E^.j. 

The linear closeness condition is then expressed as 

s < s 0 = sin(e) (8) 

If the threshold s 0 is chosen equal to zero, equation (8) defines the 
classical linear dependence. The choice of s Q greater than zero expresses 
the uncertainties that are due to real measurements and real computation, the 
latter introducing in addition some errors in the direction of the vectors 
3Y/3ai, which justifies even more strongly the necessity of the linear 
closeness concept. 

A 

In the set of the N vectors 3Y/3a^ there may be different linearly 
close subsets that will be responsible for a nonunique solution in the identi- 
fication. However, the nonuniqueness will affect only the parameters 
involved in these subsets. Therefore, if it is possible to find these 
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dependent (or almost dependent) sets, not only the cause of nonuniqueness is 
found but also some cure. Indeed, by removing from the identification process 
one parameter in each set (i.e., leaving these parameters with a constant 
value) the degeneracy of the space is removed and the identification may 

continue with a smaller number of parameters for which there will be a unique 
solution. Of course, this solution will depend upon the values given to the 
discarded parameters, but, and this is a key point, if the true value of the 
discarded parameter of a given dependent set is known from any other source of 
information, and if this value is given to this parameter, then the values 
obtained in the identification for the other parameters of this set are 
correct. Moreover, even if nothing is known, some parameters might not be 
involved at all in any dependent set, and for these independent parameters a 
correct answer will be obtained, regardless of the situation for the others. 
That is, if the dependent parameter is removed and the identification con- 
tinued, the final value of all parameters that were not involved in the 
dependent sets will be correct. 

An obvious point, which has also to be considered, is that, above all, a 
vector SY/Sa^ has to be nonzero. If such a vector is zero, this means the 
parameter a^ has no effect on the model, at least for the^particular input- 
output set considered. Again one should realize that if SY/Sa^ is small 
enough, 2 a^ will be considered an irrelevant parameter. Therefore the 
irrelevance may be checked first, while computing the vectors 3Y/3a^. 

The problem of irrelevance, which is the simplest case of nonuniqueness, 
can be related to a linear closeness situation. For example, consider a prob- 
lem with two parameters a^ and a 2 . If 3Y/3a 2 is found small enough, a 2 is 
declared irrelevant. Let us define then a new set of parameters a{ and a 2 
such that 


i » 

a i = a 2 + a± 

a 2 = a 2 a l 

it follows that 

3Y = dY__ 3Y 

3a^ 3a^ 3a 2 


3Y = 3Y_ 3Y 

3 a 2 3 a^ 3 

and the angular separation between these two vectors is 

2 The problem of defining a good criterion is delicate and it is difficult 
to give definite rules. One possible approach is to compare the relative 
change in the output to the relative change in the parameter. Also the pre- 
cision with which the computations are performed has to be taken into account, 
since the change in the last digit is not very significant. 
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s < 2| 9Y/9a 2 |/| 3Y/a ai | 


( 9 ) 


A 

which is small when 3Y/3a 2 is small. Therefore the irrelevance of a 2 is 
equivalent to a closeness between aj and a^. 

In conclusion, the analysis of the structure of the space Ejsj provides 
an accurate way of detecting the occurrence of nonuniqueness and prevents this 
occurrence from causing a complete loss since it indicates not only the 
correctly identified parameters, but also what a priori knowledge is required 
to complete the identification. How such an analysis can be performed on a 
digital computer will be discussed next. 


DEPENDENCE ANALYSIS 


Two important steps that are required before considering the computer 
program are discussed in this section. The first step is concerned with the 
definition of the operations that must be performed and the mathematical prin- 
ciples on which they are based. The second step establishes the actual method 
of computing the numerical values corresponding to these operations. The 
mathematical background used in the section "Principle" below can be found in 
appendixes B and C, while the detailed analysis of the computation technique 
is given in appendix D. The program itself, which has been written in 
FORTRAN IV, is explained and listed in appendix E with an example of output. 

Since the dependence analysis involves the N vectors 3Y/3a^, it is of 
course supposed that the components of these vectors have been computed and 
stored, either as a part of the identification process, or only for the pur- 
pose of this analysis. The actual computation of these components will not be 
discussed since the technique may differ considerably depending upon the model, 
the type of problem and other practical considerations. However, it might be 
convenient to look for the irrelevant parameters at the time of this computa- 
tion (i.e., check the magnitude of 3Y/3a^ with respect to some criterion). 

In any case, a check for zero valued vectors is always made at the beginning 
of the dependence analysis. 


Principle 

The set of the M-dimensional vectors 3Y/3a^ must be analyzed. Since 
the linear closeness involves only angular properties, a set of normalized 
vectors 


Di = (3Y/3a i )/|9Y/3a i | i = 1,N (10) 

may be used instead. At this stage the irrelevant parameters have already 
been found and the corresponding vectors 9'Y/9a. removed from the set, so 
that N is the number of relevant parameters and equation (10) is meaningful. 
Call the set of the N vectors 9Y/9ai- The problem is to find what are, 

if any, the almost dependent subsets Bj( in the set B^ (Bj( C B^, K = 1,N). 
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One possible method is the following: since an y subset that contains a 

linearly close subset is itself linearly close (appendix B) , one could first 
look for all the possible two-dimensional sets, then for the three-dimensional, 
and so on up to N-dimensional . Each time a subset is found linearly close, 
one vector could be removed and the process could be continued with a smaller 
number of vectors. This way would ensure that when a K-dimensional subset is 
found dependent (or almost dependent), by removing one vector of the subset, a 
K - 1 dimensional subset would be obtained that would contain only independent 
vectors. Unfortunately, the time involved makes the straightforward applica- 
tion of this method to a digital computer unfeasible since the computation of 
2^ determinants is required (which already takes one minute on the IBM 360/60 
for 12 parameters and this figure is doubled for each parameter added). 

Another important consideration is that when a vector has to be removed, 
many choices are possible and some are better than others. We have seen 
indeed that the nonuniqueness comes from the fact that the basis vectors 
BY/8ai constitute a skew basis. Conversely, the more nearly orthogonal the 
basis, the less tendency to nonuniqueness. For example, if three vectors are 
found in a plane (fig. 9), the best choice is to discard D 2 because it 
leaves D]_ and D 3 , which are almost orthogonal, to construct a two-dimensional 
basis. The worst choice is to discard Di which leaves a skew basis. There- 
fore it is advantageous if the technique provides also a way to optimize the 
basis. This can be done by constructing systematically an optimal basis with 
the vectors D±. To start the process a first vector is chosen. Call ni 
the index of this vector (e.g., if D 3 is chosen, m = 3) ; the parameter a^ 
corresponds to the first vector D ni • Many ways, optimal or not, may be 
devised to choose this vector, but generally it is a matter of common sense to 
know which parameter at least should be kept in the model. The next step is 
to find which of the N - 1 remaining vectors is the farthest from • To 

do that, the distances si from the Di to D ni are computed (fig. 10). If 
s n2 is the largest, then the corresponding vector D n2 is chosen to con- 
struct a two-dimensional subbasis with Dn^ • Also, the smallest value of the 
si is searched and compared to the threshold s Q . If it is larger than s 0 
the name of the corresponding parameter is kept in memory as well as the value 
of si. This defines a "critical" parameter, in the sense that, if the thresh- 
old was increased, it would have been considered dependent. Therefore the 
knowledge of this critical value si might give an indication of potential 
difficulties. If Si is less than s 0 , the corresponding parameter is 
declared dependent and Di removed from the set of vectors. The same process 
is now applied to find the third basis vector D n3 - That is, all the dis- 
tances from the N - 2 remaining vectors to the plane (Dm* D n2 ) are computed, 
and D n3 corresponds to the largest distance. On the example of figure 10, 
when the four-dimensional subbasis is being formed, the distance from Dj to 
the three-dimensional subbasis (D ni > D n 2 > D n 3 ) is found to be less than s 0 . 
Therefore Dj is discarded as linearly close to this subbasis. To memorize 
this event involving the parameter aj , a corresponding integer MN(j) is set 
to 0 in the program. This integer, or "dependency index," will be used 
internally by the computer to keep track of the dependent, independent or dis- 
carded parameters. Then, by computing the projections of Dj on the subbasis, 
Dj is found to be in the plane (D n2 , D n3 ) . Therefore the actual almost depen- 
dent set is (Dj , D n2 , D na ) . This is memorized by setting MN(n 2 ) and MN (n 3 ) 
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equal to 1 (meaning that a n2 and a n3 are dependent parameters) and by 
storing j ,n 2 and n 3 in a logical array (equivalent to the information, the 
first dependent set contains the parameters aj^ a nz and a n ^) . More details on 
the program are given in appendix E. 

The above process will continue until all the N vectors have been 

either used in the basis or discarded. Since MN was set to 2 for all the 

parameters to begin the analysis, if a vector has never been involved in any 
dependent set, its dependency index will keep this value, so that a 2 will 
indicate an independent parameter. Finally a union of the dependent sets may 
be performed if they have a nonempty intersection. For instance if the depen- 
dent sets (a 5 , a 7 ) and (a 2 , a 5 , ag) are found, they may be combined in the set 

(a 2 , a 5 , a 7 , ag) , and in this case if both a 7 and ag are discarded, the 
final values found in the identification for a 5 and a 2 will depend upon 
those of a 7 and a 8 . It is essential to note at this point that if some 
linearly close subsets have been found, the final basis will be composed of 
only N ! vectors (N ! < N) . This means that only N’ parameters can be 
uniquely identified. Hence N ! - N parameters have to be removed from the 
identification process (i.e., their value will not be changed). An optimal 
choice is to remove those for which MN = 0, and this can be done automati- 
cally in the program. Otherwise the user has to decide which parameter he 
wants to remove in each dependent subset and reset his identification program 
correspondingly. Therefore, although the removal of some vectors Di is 
required in the dependence analysis in order to find the dependent subsets, it 
does not imply the removal of the corresponding parameters. That operation 
belongs only to the application of the analysis. 


Computation Technique 

Once the irrelevant parameters have been eliminated, it is possible to 
compute the normalized vectors and start the construction of the optimal 

basis. At each step of this construction, the projections of a vector on the 
last subbasis must be computed. This still seems very strenuous, but we have 
found that all these quantities can be obtained during the recursive computa- 
tion of a single N x N determinant (which is quite an improvement compared 
to 2 n determinants needed in a straightforward method) . 

Consider indeed the matrix D whose N columns are the M-component 
vectors D-[ and form the product 

G = D t D 

where G is an N x N matrix, positive semidefinite, and the elements g-jj 
are equal to the inner products dTd.. 3 Because of the normalization, the 
diagonal elements are l's. The determinant of G is called the Gram deter- 
minant of the vectors Dj , D 2 , . . . , Dj\j, and is equal to the square of the 

3 It is also interesting to have a quick look at this matrix because if an 
off-diagonal element gjj is equal or very close to 1, it indicates that Dq 
and Dj are very close (i.e., the parameters a^ and aj are dependent). In 
this case s is simply equal to (1 - g 2 j) 1/2 . 
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volume of the N-dimensional hyperparallelepiped determined by these vectors 
(ref. 9) . 

DET(G) = 

Computing this determinant is equivalent to computing a volume. Since 
there is a recursive method of computing an N-dimensional volume which 
involves the distances s^, a recursive method of computing the determinant 
should exist that involves these quantities. If is the volume corre- 

sponding to a K-dimensional subbasis , the volume obtained by adding a unit 
vector Dj^ + i, at the distance s^+i from the subbasis, is indeed 

V K+1 = V K x S K+1 


Squaring the two sides of this relation leads to the recursive relation for 
the determinant. It is shown i-n appendix D that this can be obtained by 
pivoting about the diagonal elements. At the K th step, the situation is as 
represented in figure 11. By rows and columns exchange, and renormalization, 
the K first rows (corresponding to the K basis vectors D ni , D n2 , . . ., 
Dp^) form an upper triangle with the diagonal elements equal to 1. The other 
diagonal elements such as gj^, gjj, • • • are, respectively, equal to the 
square of the distance from , Dj , . . . to the subbasis (D ni , D n2 , . . ., 
D n j() • Indeed, it can be seen directly that the determinant corresponding to 
this subbasis plus , say Di, is equal to the determinant corresponding to the 
subbasis multiplied by gii* Symbolically 

det K+ 1 = det K x gii 

which is precisely the recursive relation discussed above. Therefore search- 
ing the elements gii> gjj> * • • f° r the largest and the smallest will end 

in the determination of the next basis vector and of the critical one. If for 
instance gjj is found to be less than s q , then it is a simple matter to 
obtain the components of Dj on the subbasis to determine the members of the 
dependent set (appendix D) , because the system of equations to be solved is 
already triangular. In the same way a few extra computations at the end of 
the analysis gives the solution dA f of the equation 

(D t D)dA' = D T ER 

from which the solution for equation (5) is simply obtained by 
da t = -daplaY/SajJ (i = 1,N and MN (i) j- 0) 


Computer Output 

During the computation of the Gram determinant the following information 
was obtained: 

1. Name and separation of the new basis parameter 

2. Name and separation of the critical parameter 
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3. Names of the parameters found in a dependent set 

4. Name and separation of the parameter optimally discarded from this 

set 

5. Also correlation coefficients between this parameter and the other 
parameters of the dependent set (see appendix A) 

6. Values of the dependency index (0, 1, or 2) for each parameter 

7. United sets and name of the optimally discarded parameters 

These are the type of items an engineer may want and this is what we have 
displayed in our program (appendix E) . 

Certainly, there are as many ways to use this subroutine as there are 
users, but some hints may still be indicated in this discussion* Two kinds of 
corrective action can be undertaken following this analysis. One is entirely 
automatic and controlled by the zeroes of the dependency index that are used as 
a signal for the computer to bypass any computation related to the correspond- 
ing parameter. This is compulsory if a matrix inversion is used in the iden- 
tification algorithm. It will speed up the computation in other cases. The 
second kind of action is a human decision based on the information obtained 
and a general insight of the particular problem. Consider for instance the 
list of the basic parameters. It will in fact contain all the parameters if 
the subroutine is called with a threshold equal to 0.0. Because of the opti- 
mization procedure, the parameters are sorted by increasing degree of depen- 
dence. Since the separation is also a measure of the confidence in the 
parameter estimate (appendix B) , this list expresses the intrinsic properties 
of the identification independently of any statistical property of the noise. 
Consider now the choice of the discarded parameter when a threshold is speci- 
fied. What is indicated by the computer is only an optimal choice, which will 
give the best results in the identification of the remaining parameters. In 
the example in appendix E, a g and a 14 were discarded in the first dependent 
subset. If the value of these parameters could be obtained from other measure- 
ments, the identification could be run again with these values. However, the 
engineer may find that known values are available only for ag and a^, in 
which case he will decide to identify a g , an, and ai 4 . This kind of trade- 
off is typical of real problems and cannot be ignored. It requires more skill 
than brute logic; therefore it cannot always be mechanized automatically on 
the computer. 


APPLICATION TO THE IDENTIFICATION OF A SATELLITE SIMULATOR 


Data Acquisition System 

The Ames Research Center three-degrees-of-freedom satellite simulator is 
an air-bearing supported platform provided with a system of servo-gimbals and 
mercury connectors allowing accurate measurements of the attitude angles and a 
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torque-free link with the outside (ref. 10). This particular air-bearing 
configuration allows about ±20° amplitude in roll and pitch and there is no 
limit in yaw. 

On figure 12 is sketched the general configuration of the whole system 
for data acquisition and processing used for the identification of the plat- 
form. The attitude rates (p, q, and r) are sensed by three rate gyros whose 
analog outputs are digitized, after filtering, and recorded on magnetic tape 
by an E.A.I. 8400 computer. This tape has to be converted to an IBM com- 
patible tape from which data cards can be punched. These cards contain there- 
fore the dynamic data (i.e., the discrete time histories of the three rates). 
The initial attitude angles of the platform (Euler angles) are read on the 
gimbal readout display and the corresponding values punched manually on cards 
(static data). The identification is run on an IBM 360 digital computer which 
receives the identification program. A special set of control cards is used 
to specify the different operations performed by the identification program. 


Equations of Motion 


The nonlinear equations used to describe the system in the case of large 
angular motion may be written: 


g = S (aj)g 




(11) 

h = S(to)h + T 




(12) 
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where ca is the angular velocity of the platform (of components p, q, and r) , 
S(w) is a matrix expression for the crossproduct: 


( 0 r 

r 0 

q -p 



(S(u>)g = g x co) 


g is the gravity acceleration vector, h the angular momentum, J the inertia 
tensor of the platform and T represents the applied torque. 


These equations have been developed in reference 11 for attitude control 
analysis. They are written here in body coordinates with the origin at the 
center of the bearing. The first equation is purely kinematic, as it 
expresses the change in the coordinate systems (in its general form it 
involves the whole attitude matrix, but in this particular problem, where no 
control in yaw angle was considered, the direction of the vertical only was 
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required), the second equation is a torque equation, and the last one is just 
the classical definition of the angular momentum. 

In normal operation (i.e., for satellite simulation purpose) the platform 
is completely balanced, that is, the center of gravity coincides with the 
bearing center. In our experiments the center of gravity was offset deliber- 
ately to make the system pendulus. If M is the total mass and b/M the 
bearing offset (three components vector), the torque exerted by the gravity is 

T G = -S(b)g (14) 

The natural damping of the system, mainly due to aerodynamic effects, was 
quite small and in some tests it was convenient to use the reaction wheels as 
dampers. The three axes were damped separately; the corresponding control 
system is shown on figure 13 for one axis. Considering one axis, we use the 
classical linear model for the DC motor, 

E = RI + KgQ (15) 

where E is the applied voltage, I the current in the motor, R the resis- 
tance of the wiring Kg a back electromotive force constant, and Q the 
angular velocity of the shaft (which is also that of the wheel). The torque 
available is proportional to the current and satisfies 

T M = Kpl = Jwft + Kq5> (16) 

where is the inertia of the rotor-wheel system and Kp a damping con- 

stant. When the wheel is used as damper, the applied voltage E is con- 
trolled by the platform rate, say the roll rate p: 

E = C lP (17) 

where C^ is a control gain. Eliminating E and I from equations (15), 

(16), and (17) and calling h^i the angular momentum of the roll wheel 
(h W i = J^Q) , we obtain finally an equation of the type 


h Wl " c lP " a l h Wl 


(18) 


where is a constant. 

Considering now the total angular momentum h^ of the wheels (of com- 
ponents hwi , hw 2 , hw 3 ) 5 the control equation may be written as 

h w = Cm - ah w (19) 

where C and a are diagonal matrices corresponding, respectively, to control 
gains and back electromotive force plus damping factors. The main nonlinear- 
ities were introduced by the torque and speed limiters. They enter the 
preceding equation as limits in the values of the components of Cm and h^, 
respectively. 
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Finally the torque exerted by the wheels on the platform can be expressed 
as 


Tjv - s (o))bw - h w (20) 

and the complete system is modeled by letting the torque T in equation (12) 
be 


T = T w + T g (21) 

When the wheels are not used, T \\j is replaced by a linear term in w to model 
the natural damping, that is DMPuj, where DMP is a diagonal matrix. 


In all cases, the "output" of the system was obtained by the three rates 
p, q, and r seen through the first-order filters and the digitizer (fig. 12). 
Thus three equations have to be added to model the system completely: 


xyi = k x p - y 1 ' 
xy 2 = k 2 q - y 2 . 
xy 3 = k 3 r - y 3 , 


( 22 ) 


The DC gains of the filters (k) were obtained by independent calibration 
of the gyros and computer links, and the time constant t was known from the 
setting of the analog filters. 


Finally the complete sets of equations corresponding to the free pendulus 
(natural damping) case, and the reaction-wheel-damped case may be written: 


Case A: natural damping 
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with the equations of the filters: 


= kip - y 1 
xy 2 = k 2 q - y 2 
T y 3 = k 3 r - y 3 

The parameter a 13 appearing in the last torque equation was needed to 
describe a constant torque in yaw experienced during the tests (because of 
some drift of the system not corrected at this time). 

Case B: reaction-wheel-damped case 




Ty x = k x p - y 2 

xy 2 = h 2 q - y 2 

xy 3 = k 3 r - y 3 


with the limiting conditions: 

(motor torque) | C x p | < a 16 , |C 2 q| < a 17 , |C 3 r| < a 18 

(wheels spin) |h Wl | < a 22 , |h W2 | < a 23 , |h W3 | < a 24 
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Parameter Identification Problems 


An analysis of the equations can show that, in case A or B , there is not 
a unique solution for all unknown parameters. One parameter has to be known; 
hence the lateral offset b x was measured independently because it was easy 
to obtain by adding a small weight on the side of the platform, at a known 
distance from the bearing. For case A, 12 unknown parameters remain to be 
identified: the 6 inertias and products (Jx, Jy, Jz, Jxy, Jxz, Jyz) , the 

2 center-of-gravity offsets (by, b z ), the 3 damping terms (DMP X , DMP y , DMP Z ) 
and the turbine torque a^ 3 . In case B, 21 parameters are to be identified: 

9 of the parameters as in case A (all except DMP X , DMPy , DMP Z ) plus the 3 con- 
trol gains Ci, C 2 , and C3, the 3 back electromotive force and motor damping 
terms 04, a 2 , a 3 , the 3 limits in motor torque a^g, a i7> a l8> and the 3 
limits in wheel momentum a 22 , a 2 3, a 2 4. 


Experimental Results 

In the two cases, the platform was first held in some position by an 
electromechanical device, and the attitude angles (pitch 0 O and roll tj> 0 ) 
measured. The initial conditions therefore were 

g 1 = -sin(0 o ) 

g 2 = cos (0 o )sin(cf) o ) p = q = r = 0 
g 3 = cos (0 o )cos (cj) 0 ) 

Because of the gravity restoring torque, once the platform was released 
an oscillatory motion was observed about the three axes, and the corresponding 
discrete time histories of p, q, and r were recorded. About 200 points were 
used in each, covering 5 or 6 periods, thus the vector Y had about 600 
components . 

When the motion of a body is excited about one axis only, it only depends 
upon the inertia about this axis and no information is available concerning 
the other inertias. This has actually been the current method to determine 
the inertias of aircraft where great care is taken to obtain single-axis exci- 
tations. In the case of the platform, it is intuitive that motions about the 
three axes should be excited to obtain a good identification of the inertias 
and products. Therefore some difficulties may be expected because the yaw 
motion is poorly excited by the gravity restoring torque. This motion depends 
of course on the lateral center-of-gravity offset and the initial conditions, 
but both are limited by the angular limitations of the platform, and in our 
tests the yaw rate was at most 0.1 or 0.2 of the roll or pitch rate. 

Case A- These data were taken in the free-oscillation case (no reaction 
wheels damping), with the initial conditions: 

p, q, and r = 0 
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ii min mi 


pitch angle = -8.44° 
roll angle = 10.44° 

Lateral bearing offset on the x axis: b x = 0.211 ft- lb. The first identi- 

fication was run with a threshold set at s 0 = 0. All 12 parameters were 
identified and a very good fit in the time histories was obtained after a few 
iterations. This fit was measured by a quantity denoted as "COST" which was 
simply the sum of the square of the errors, that 'is |ER| 2 . The decimal loga- 
rithm of this quantity is plotted on figure 14 as is the evolution of the 
values for the six elements of the inertia tensor. Although the behavior of 
the parameters is somewhat erratic at the beginning, they finally reach an 
asymptotic value for which the cost remains essentially constant. Therefore 
one might have been satisfied with these results and considered that these 
final values of the parameters were correct. The dependence analysis, however, 
indicated a low value for the separation of the parameter (J xz ) . Consequently, 
the identification was run again, with the same initial values for the param- 
eters, but with a threshold set at s 0 = 0.02. As it appears on figure 14, 
the evolution of the parameter values is quite different in this run. Indeed, 
J xz was discarded at the second and following iterations because of its depen- 
dence on J x , J z , and ax 3 . The diagnostic issued by the computer when the 
final values were obtained is given in appendix E. The important fact to note 
here is that these final values are different from those obtained before, yet 
the match in the time histories is as good. To illustrate this similar match, 
a plot of computed and measured time histories of the roll rate is shown on 
figure 15 for these two identification runs. 

The final numerical results obtained in these two runs are given in 
table 1. These results deserve some comments and explanations. The "error 
bounds" were computed according to equation (B 8 ) : 

( da K > 2 " 2, ,2 

s K |3Y/3a K | z 

where the measurement error e 0 was estimated from the residual error between 
model and plant output when the best match was obtained. Note that the values 
of the independent parameters (dependency index = 2 ), although different in 
each run, are within the predicted error. The dependent parameters, however, 
(dependency index = 1 ) exhibit large variations as a result of the nonunique- 
ness. Another quantity given in these tables is the sensitivity . It is a 
measure of the model response to parameter change and it is defined by 


Sensitivity to the parameter ai 



3Y 
3 a-^ 


A sensitivity of 1 corresponds to a direct proportionality between the corre- 
sponding parameter and the output. If the sensitivity is very small with 
respect to 1 it indicates that the output is almost independent of the param- 
eter, and large errors may be expected in the estimation of this parameter. 
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Therefore this quantity can be used as a criterion for deciding whether a 
parameter is irrelevant. Sensitivity and error bounds are not a part of the 
dependence analysis but are computed independently in the identification 
program. 

Finally, a partial check of these identifications was made. The value of 
b z , vertical offset of the center of gravity, was determined independently by 
a static test of the platform, and found equal to 0.78 ft- lb, which agrees 
well with the results of the two identification runs. 

Case B- To test the algorithm in a strongly nonlinear case, another 
identification was performed with reaction wheel damping. It provided an 
interesting illustration of the problems encountered and how they were 
detected by the algorithm. 

The initial conditions for this case were: 

p, q, and r = 0 
pitch angle = -0.65° 
roll angle = 12.98° 

Lateral bearing offset on the x axis: 

b x = 0.134 ft- lb 
threshold = 0.01 

On figure 16 are shown the computed and measured time histories of the roll 
rate. (The wheel spin is also plotted exhibiting sharp saturation effects.) 

A residual systematic (not random) error is observed here that must be attri- 
buted to a modeling error. Indeed, many assumptions were made in modeling the 
control system, mainly concerning the linearity of the DC motors and the 
absence of friction in the wheel bearings. However, considering the approxi- 
mations made, the agreement is quite satisfactory. The final results of the 
identification are given in table 2. The parameters a 18 and a 2 4 were found 
to be irrelevant. The latter represents the limit of the yaw wheel momentum, 
corresponding to the spin limitation of the yaw wheel. Indeed, the yaw wheel 
momentum time history showed that it never exceeded its limit; therefore it 
was not possible for the computer to determine this limit. For a similar 
reason, the limit in torque a 18 could not be determined. (Note that the 

sensitivity of these two parameters was found equal to zero.) 


CONCLUDING REMARKS 


A method has been described that determines which parameters of a system 
model can be uniquely identified from a given set of measurements. It applies 
to nonlinear systems and does not require any special form for the model. 
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providing that the model is completely defined and depends upon a finite 
number of parameters. The basis of this method is that all the necessary 
information is contained in the derivatives of the error vector with respect 
to the parameters. These derivatives are vectors in an M-dimensional space 
defined by the measurements (their components are often called "sensitivity 
functions" or "influence functions"). The nonuniqueness has been shown to be 
the result of linear dependence between these vectors in the ideal case of 
perfect measurements and exact computation. The concept of "linear closeness" 
has been introduced to take into account the uncertainties introduced by the 
errors occurring in the measurements and in the computation, thus allowing a 
realistic analysis of practical problems. 

To perform this analysis, an algorithm has been developed for digital 
computation that can easily be integrated in a complete identification proce- 
dure. This is done practically at no extra cost in computation time if the 
procedure already requires the computation of the sensitivity functions (as in 
multiple linear regression or in quasi linearization for instance). 

The detailed comments put out by the computer not only point out where 
the identification has, or is going to fail, but also indicate what could be 
done to remove the nonuniqueness. A very useful feature of this algorithm is 
its ability to obtain the correct values of some parameters (independent param- 
eters) despite the fact that the complete solution might not have been unique. 

Successful results obtained in the identification of a three-degrees-of- 
freedom satellite simulator indicate that the air-bearing technique might be 
applied to the determination of the inertia tensor of an actual aircraft. The 
identification of a highly nonlinear control system implemented in the simula- 
tor has also proven the effectiveness of the algorithm in such cases. 

It must be noted finally that beyond the strict application of this 
technique to identification, there are interesting potentialities in other 
fields, in modeling for instance. Indeed, the algorithm indicates how the 
model could be simplified since the most dependent parameters can be elimi- 
nated by raising the value of the threshold. As the threshold is raised more 
and more, the model becomes more and more simple, although less and less accu- 
rate. Another possible use is the practical tailoring of the input to obtain 
optimal results in the identification. In this case, one could maximize the 
separation of a given parameter, or some average upon a group of parameters, 
by adjusting the parameters defining the input sequence. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, May 14, 1971 
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APPENDIX A 


LINEAR DEPENDENCE AND FUNCTIONAL RELATIONSHIPS 


The linear dependence between the K (M-component) vectors 
9Y/9a 2 , . . . , 9Y/9ay may be expressed by 


K-l 

9Y_ 9Y 

9aj( ' ■> 9a^ 

i=l 


9Y/9a! , 


(Al) 


where the Xi are some nonzero functions of the parameters a and the K-l 
vectors 9Y/3a^ are supposed independent. When a^ are varying, Y is con- 
strained to stay on a K-l dimensional hypersurface. It is generally 
possible to find a system of coordinates (zj, z 2 , • . Zk- 1^ on ^is surface 

so that Y may be expressed as a function Z of the coordinates z: 

?(ai, a 2 , . . ., a K ) = Z(zi, z 2 , . . ., z K-1 ) (A2) 

A 

Therefore Y is not really a function of the K parameters but depends 
upon K-l functions of them . If the partials of both sides of equa- 

tion (A2) are taken with respect to a and compared to equation (Al) , it can 
be shown that the functions z satisfy the system of equations: 



9a K 




3=1, K-l 


(A3) 


Given the values of a, the values of z are well defined. The converse is 
not true, but if a K is known, then there must be a one-to-one correspon- 
dence between z and the first K-l parameters, since Y is now a func- 
tion of these K-l parameters only. The necessary and sufficient condition 
on the functions z is that their Jacobian with respect to the K-l a is 
not zero; that is. 


9 (Z], , z 2 , . . ., zk-i) Q 

3(ai, a 2 , . . ., aK-i) 


(A4) 


We now examine the case where the vectors 3Y/9ai are only linearly 
close (see appendix B for definitions) . Equation (Al) becomes 


K-l 

3Y_ _ , 9Y 

3a K " 2^ 9^f + EK 

i= 1 


(A5) 
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where is a (M- component) vector orthogonal to the space Ej(_i spanned by 

the K - 1 vectors 9?/9ai and of small magnitude. In this case Y depends 
also upon a^ and may be written as 

A 

Y(ai, a 2 , . . . , aj() = Z(zi, z 2 , • • • , zr-i, a K ) 

if z± are chosen in such a way that they verify equation (A3) , then using 
equations (A3) and (A5) one obtains 


9Z 

9^ ~ £ k 


(A6) 


A 

The change in Y when a^ 


is changed while zi are kept constant is thus 


dY = dZ = £j(da^ 


A 

For a small enough e k, dY will be negligible and, because of equation (A6) , 
the property equation (A2) will be verified. 

Let us now evaluate the total change in Y for arbitrary da, 


A 

dY 



da K 


Eliminating 9Y/9a K with equation (A5) it follows that 

K-i 

dY = (dai + Xida K )r|~J + e K da K ( A7 ) 

i= 1 

If e K is small enough, then for a given change da^ there exist da^ such 
that d¥ stays close to zero. From equation (A7) we see that the values of 
such da^ are gi ven by 


dai — -Xida^ 


(A8) 


Therefore if a value was given to a^ with an error da^, the corre- 
sponding errors in the values of the ai obtained in the identification are 
given by equation (A8) . Consequently, it is very useful to compute the quan- 
tities Ai . In the algorithm (see appendix D) a normalized set of vectors is 
used : 
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OY/aaO 

D = J_ j = 1, K 

| 3Y/3aj | 

and the quantities y. are computed. The y^ are defined by transforming 
equation (Al) into 

K-l 

°K = S ^i D i 
i=l 

They are the components of Dj( upon the vectors and are related to Ai 

by 


Uj 1 9Y/3a K | 
| 3Y/3ail 
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APPENDIX B 


DEFINITION AND PROPERTIES OF LINEAR CLOSENESS 

Consider a set By of K unit vectors (Di , D 2 , . . . , Dy) spanning a 
space Ey. The first K - 1 vectors span a space EfC-i- Let ey (ey £ Ey) 
be a unit vector orthogonal to Ej^_ ^ . Then Dy may be decomposed as 


K-l 

D K = s K e K + 2 Pj D i (Bl) 

i=l 

The Pi are the components of Dy on the D^ in the basis (ey, Di, D 2 , 

. . . , Dy_j) and Sy is a positive scalar (this can always be done by a 
proper choice of ey) equal to the sine of the angle between Dy and its 
projection on Ey.^. When sk is zero, equation (Bl) defines a strict linear 
dependence between the vectors D. If sy is arbitrarily small, these vec- 
tors become almost dependent and sy is a measure of their degree of depen- 
dence or of the separation between Dy and the space Ey_i> We may therefore 
define the linear closeness as: 

The vectors of a given set B% are linearly close with respect to 
the threshold s Q iff there exists a vector Dx of the set such 
that its separation sx from the others is less than or equal to 

Sq. 

We may write this symbolically: By is LC/s 0 iff 3 Dy € By such that 

s K < s 0 . 

What has been defined is obviously an angular property, therefore it 
will not change if each vector D is multiplied by a different scalar. For 
any set of vectors, it is possible to define a corresponding set of normal- 
ized vectors. By, and iff By is LC the original set is also LC. For 
convenience and ease in the proofs, we always work with the normalized set. 

Consider a set By where Dy is given by equation (Bl) . Introduce 
another vector, say Dy+i, to construct a new set By +1 (By +1 = By U Dy + i). 

Then Dy may be written 


D 


K " 


K- 1 . 

s K+l e K+l + ? y i D i 
i= 1 L 


y K+l°K+l 


(B2) 


where ey +1 is a unit vector orthogonal to D^ and Dy +1 , and are some 

scalars. Recalling that 
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T _ , 
e K e K “ 1 


e K°l ' ° 


il'K.l " 1 e Ll D i * 0 


i - 1, K - 1 


e K« D K + l - 0 


one obtains from equation (B2) 


e K + l D K = s 


and from equation (Bl) 


therefore 


K+l 


eJ +1 D K = SKCe£ + 1 e K ) 


S K+1 = (e K e K+l )s K 


( e ^ e K +i ) is positive (because s^ and s^ +1 are positive) and less or equal to 

1 (because equal to cos(0), where 0 is the angle between the two vectors 
e) , thus we obtain the important inequality 


S K+1 < S K 


(B3) 


Property 1 - This result shows that the increase in dimensionality 
always increases the dependence (unless the new vector added is 
orthogonal to the others in which case nothing is changed) , which 
means that when the nwriber of parameters is increased there is 
more chance for a nonunique solution . 

Now suppose that a set B^ is given in which the subset B^ (B K C Bjq) 
is LC. From the definition 


SK < s 0 

Starting from B^, add the vector to construct the set Bj^ +1 C B^. 

Because of equation (B3) 

S K+1 ^ s o 

therefore B^ +1 is also LC. Continuing to add vectors this way, we see 
that ultimately sjyj ^ s Q and we have 

Property 2- If a set includes a LC subset 3 it is itself LC 3 and 
its separation is less than or equal to that of the subset. 


29 



I II I III 1 1 III 


To summarize the preceding results, let us introduce this notation: if 

B is a set of vectors spanning some space E, and D a vector not belonging 
to the set B (but belonging eventually to the space E) , the separation s 
between D and B is written as 


s = S (D/B) 


(As has been seen, 0 < s < 1.) 


Definition of the Linear Closeness 

B is LC/s 0 iff 3 D 6 B such that S(D/B') ^ s 0 , where B = B' U D, 
Properties - 

LC 1 : 


S (D/B UD')< S (D/B) V B,D and D' 

LC 2: If B' C B and B' is LC/s 0 , then B is LC/Sq with s' Q < s Q . 

LC 3: 


S (Dj/B U P 2 ) _ S (D^B) 

S(D 2 /B u Dj] S(D 2 /B) 

The proof of this last property is very simple; consider the volume of 
the hyperparallelepiped constructed with the set B U D 2 U d 2 : 

V (B UDj U D 2 ) = V (B U Di) x S(D 2 /B U D x ) = V(B U D 2 ) x SfDi/B U D 2 ) 


In the same way 


V (B U D x ) = V (B) x S (D ! / B) and V(B U D 2 ) = V(B) x S(D 2 /B) 


which leads to LC 3 by substitution. 

Relationship Between Separation and Error in Parameter Estimates 

Let s k be the separation between D^ and the set Bj(_i. We have as 
in equation (Bl) 

K- 1 

°K = s K e K + 12 Pi D i (B4) 

i=l 
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When there is some ,! noise n in the measurement vector Y, not only Y 
lies outside of the subspace E^, but also its projection Y Q on this sub- 
space is shifted by some amount SY 0 from its value in the absence of noise. 
As a matter of fact, because the identification procedure tends to make the 
computed vector Y reach the target Y 0 , the component of the noise orthog- 
onal to E k does not affect the result, whereas the component in Efc, SY 0 , is 
directly responsible for the error in the parameter estimates. This uncer- 
tainty 6Y Q may result from random perturbations of the system, or instrument 
errors, quantization errors, nonlinearities or other unknowns, but we may 
assume that some bound is known for 6Y 0 , for instance. 


< • s *° < ^ 

Therefore a set of parameter values that causes Y to satisfy 
dY T dY = (Y - Y 0 ) T * (Y - Y q ) < e* 


(B5) 


(B6) 


is as good (or as bad) as any since Y 0 is not known exactly. If the param- 
eters are changed by the amounts da^ , da 2 , . . da^, the corresponding 
change in Y may be written as 


dY 


K- 1 


i=l 


A 

3Y 

D-; 

da^ + 

3Y 

3a^ 

1 


3a. K 


°K da K 


using equation (B4) leads to 
K-i 


dY 


i=l 


8Y 


3Y 

, \ n 

3Y 

3a^ 

da i + H 

3a K 

da K ) D i + S K 

8a K 


da K e K 


(B7) 


or 


dY = V + ( s , 


( S K :fi£ da k) e K (where VG E K-1 and e K _L V ) 


finally 


dY T • dY = |VI“ + s 


2 . „2 


3Y 


8a K 


(d a K ) 2 


Comparing with equation (B6) gives us an upper bound for da^, the error in 


(da K ) 


2 < 


2 

e O 


S J|3Y/3 a K r 


(B8) 


31 



When ay is not coupled with the other parameters (i.e., when Dy is 
orthogonal to Ey_j), Sy is equal to 1. Thus the quantity e^/|3Y/8ay| 2 
may be interpreted as an "uncoupled error." The effect of the closeness is to 
increase the error in the parameter estimate as appears clearly in relation 
(B8). 


It is interesting to compare expression (B8) to the classical least 
squares result for the variance of the parameters (ref. 2), obtained also in 
quasi- linearization methods (ref. 12): 

E(6a|) = (DER T DER)**E(6Y 2 ) (B9) 

(DER^DER) * is the k th diagonal element of the matrix (DER^DER) 1 and it 

KK ? \ * I 2 

can be shown easily that it is just equal to l/(s£| 9Y/3ay| ). But, whereas 

(B7) supposes a white ergodic process for 6Y, our result assumes only the 
existence of some kind of bound on the measurement error. 
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APPENDIX C 


DEFINITION AND PROPERTIES OF THE OPTIMAL BASIS 


Given a set of vectors belonging to a K-dimensional space E K 
vectors of this set Dj, D 2 , . . . , D K are said to form an optimal basis in 
E{( with respect to the threshold s 0 if they have the following properties: 

For any i (i from 1 to K - 1) and any j (i + 1 ^ j < K) : 


OBI 

SCDj/Bi) > s 0 

0B2 

S(D i+ i/Bi) > SCDj/Bi) 

where B^ is defined by 



Bf = D l U D 2 U D 3 . . . U Di 

Given a set of N vectors (D]^, D 2 , . . ., D^) and a threshold s 0 it 
is possible to construct an optimal basis using the recursive process: 

(a) To start the process, an arbitrary vector, Dj, is chosen from the 

set. 


(b) The following steps are then taken: 

Step 1 B]_ = Di 

D 2 is the vector of the set such that S (D 2 /B 1 ) (Di/B]J 

Vi t 1,2 


Step j Bj = Bj _ x U D j 

Dj +1 is the vector of the set such that S(Dj + i/Bj) ^S(Di/Bj) 

Vi f 1, 2, . . ., j + 1 


(c) At each step j , D^ is discarded from the set if S(D^/Bj) 

It is clear that (c) leads to the property OBI and (b) to 0B2. 
some vectors might have been discarded, K is less than or equal to 


Since 

N. 
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We shall now derive some properties of the optimal basis from its 
definition and the properties of the linear closeness shown in appendix B. 

0B3- In an optimal basis, the vectors are sorted by decreasing values of their 
separation from the preceding subbasis: 

S(D i+ i/Bi) < SCDi/Bi_!) 

Proof : 

S(D i+1 /Bi) = SCDi+i/Bi.! UDi) < SCDi+i/Bi.j) (from LC1) 

and 

S CD i+1 /Bi_ i) < S (Di/Bi. j) (from 0B2) 

0B4- In any subbasis Bj 

S(D j _i/B j _ 2 U Dj ) > S(Dj/Bj_ 1 ) 

Proof: From LC3 and S(Di/Bj_ 2 U D j-i) = S(Di/Bj__ 1 ) 

s ( D i_l/Bi- 2 ) 

S(D i _ 1 /B i . 2 UDj) = S(D j /B i _,) x J J 

3 J 3 S (Dj/Bj _ 2 ) 

The result follows because the fraction is larger than 1 because of 0B2. 
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APPENDIX D 


COMPUTATION TECHNIQUE USED IN THE DEPENDENCE ANALYSIS 


This appendix describes the elementary operations performed by the 
computer when executing the dependence analysis program whose FORTRAN IV 
version is given in appendix E. 

Given the N (M- components) vectors 9Y/3a^, first the normalized vectors 
are computed 


Dj. = (3Y/3a i )/|3Y/3a i | 

Then the elements gij of the Gram matrix G corresponding to the N 
vectors Di are computed 


(All the diagonal elements are unity and the determinant of this matrix is 
DET(G).) Consider now the operations: 

1. Define the quantity 


DETi = g n = 1 


Multiply row 1 of the determinant successively by the elements of 
column 1 and subtract from the corresponding rows to obtain zeroes in column 1. 
The Gram determinant becomes: 


with 


DET(G) = 


T 



12 

§13 • 

' * g lN 

f 

f 

f 

22 

§23 

g 2N 

t 

1 

! 

32 

§33 

g 3N 


g N2 g N3 g NN 

g . g . i , i = 2 ,N 

6 il 5 lj ,J 
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Now divide row 2 by g' 22 which is classically called the "pivot." Call 
it p 2 . 

2. Define 

DET 2 = p 2 x DET i = g' 22 

The determinant is now 

1 g 12 § 13 • • * Sum 

o i g 23 • • • g 2N 

1 f ! 

DET (G) = p 2 x 0 g 32 g 33 * ' ‘ g 3N 

0 g N2 g N3 ' ’ ' g NN 
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DET (G) = 1 x p 2 x P 3 


'K-l 


1 x 

0 1 

0 0 

0 0 


x x 

X X 

1 X 

X X 


X 

X 


X 

X 


/ K rows 


0 0 ... x x ... x 


where x is used just to indicate the location of nonzero elements. 

Obviously, if there were only K rows and K columns (corresponding to 
K vectors DjJ , DET(G) would just equal DETj(. Therefore DET^ in general 
is the value of the Gram determinant of the first i vectors. Since the Gram 
determinant is equal to the square of the volume constructed with the vectors, 
we may write 


DET K = V(B K )2 = VCBk .!) 2 X SCDk/Bk.-l ) 2 


DETk.! = V (B k _ ! ) 2 
Compared with equation (Dl) 


P K = SCDk/Bk.O 2 

This shows that the pivotal element is precisely the square of the 
separation of Dj( from the preceding subbasis B^_ 1 . 

The construction of an optimal basis requires that, at each step, the 
largest diagonal element be chosen as pivot. In general the sequence of vec- 
tors will not be 1, 2, 3, . . . , N as described before, but it is always 
possible at each step, to exchange rows and columns to build up the triangular 
form. In practice, however, such an exchange is time consuming and will not 
really be performed by the computer. It is replaced with advantage by a 
bypassing process. Indeed, once a diagonal element has been chosen as pivot, 
say p^, and the corresponding row normalized, row and column i will not be 
changed any more. To memorize this fact that the elements of row and column 
i must not be modified, a logical variable associated to the index i is set. 
The general formula, by which the new value of the elements g. . is obtained, 
may be expressed in the general form 1 * J 

^ g ij^new ” ^ g ij^old g i next X g next j 
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where "next” is the index of the pivotal element. Each time this formula is 
to be applied, the indices i and j are checked and this computation is by- 
passed if row i or j has been already normalized. Also, when the next 
pivotal element is searched, all the pivot values are examined from 1 to N 
except those already used since the same pivotal row cannot be used twice. 

In the same way, whenever a row, of index i say, is discarded 
(corresponding to the parameter a-jj , it is memorized and all computations 
involving this index are bypassed. This is also done, right from the beginning, 
for the irrelevant or unused parameters. 

As we have seen previously, if a diagonal element, say g^, is found 
smaller than the square of the threshold s G , the corresponding vector 
has to be discarded; that is, row and column K are not considered any more 
in the computation process. In this case a special procedure is started to 
determine the subset of B^-i to which Dj( is linearly close. This is done 
by computing the components \i± of Dj( upon the K - 1 vectors belonging to 
The following equation has to be solved for the y^: 


Ei^iDi = D K i = 1, K - 1 


Multiplying both sides by Dj gives the equivalent system: 


^i y i g ji = g jK i and j = 1, K - 1 


It can be shown that an equivalent system is obtained by transforming g. . 
with the same rules used to reduce the determinant. Since the gpj , 1 ‘* 

at the step K, form an upper triangle, it is then a simple matter of backward 
elimination to obtain the unknown up. 

In the same way, the solution of any system of the form 

D.g.. da; = Cj (i and j = 1, N and MN(i), MN(j) i 0) 


where the daj are the unknowns and Cj are given quantities, can be solved 
be incorporating the cj as an extra column to the Gram matrix and performing 
the corresponding transformations. In this case, however, the pivot will 
never be chosen in this extra column, and the backward elimination is per- 
formed once the Gram determinant is completely transformed in its triangular 
form. If cj are precisely chosen as the components of the vector D T ER, 
the daj are the components of the vector dA f solution of the matrix 
equation 

(dT D)dA' = D T ER 


38 



It is then easy to obtain the solution dA of equation (5), 

(DER T DER)dA = -DER T ER 


which is simply 

dai = -daJ/jBY/aail (i = 1,N and MN(i) f 0) 
dai =0 if MN(i) = 0 
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APPENDIX E 


A FORTRAN IV PROGRAM FOR THE DEPENDENCE ANALYSIS 


This program has been written as a subroutine activated by the statement 

CALL COR(Nl, N3, THR, MODE) 

N1 number of parameters of the model 

N3 total number of output measurements 

THR value of the threshold 

MODE integer controlling the output; may take values from 0 to 6 (0 prints 
nothing, 6 prints out the full information, between 0 and 6 partial 
results are given). If 10 is added (i.e., from 10 to 16) the inversion 
of the Gram matrix is also performed. 

The other variables needed as input or output for this subroutine are 
passed via labeled COMMON'S which are used throughout the whole identification 
program. 


1 < N1 < 36 
1 < N3 < 2400 
0 < THR < 1 




Input Variables 


MNS 

controls 
MNS (I) = 
MNS (N1 + 
equal to 

the parameters to be actually used in the identification: 

0 if the user does not wish to identify the parameter aj . 

1) determines the first vector to start the analysis. If set 
0, the subroutine will choose the first vector. 

DER 

contains 

the vectors 3Y/3a (DER(j,i) = 3Yj/3aj_) 


ER 

contains 

the error vector (ER(j) = Yj - Yj) 


MN 

is used as input, only to indicate the irrelevance 

of a parameter 


(MN(i) = 0) if this has been found during the computation of the DER's 

Output Variables (Returned to the Calling Program) 

C contains the elements of the normalized Gram matrix, optionally will 

contain the elements of the inverse if MODE ^ 10 

PIV contains the values of the pivots 

AD contains the norms of the vectors 3Y/3a 

DDA solution vector of the equation (DER^ DER)DDA = DER^ ER 
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MN magic number used to control the identification (dependency index) : 

MN(i) = 0 means a^ is dependent (or irrelevant) and will not be 
identified 

MN(i) = 1 means a^ is dependent but will be identified 
MN(i) = 2 means ai is independent and will be identified 


Output Variables Used in the Print-Out Only 

These variables are not returned to the calling program but are printed 
out at different steps of the computation. Therefore some may have different 
meanings in different parts of the program. The meanings given below corre- 
spond to those they have in WRITE statements, in the order they appear in the 

program. 

NEXT index of the next basic vector found to build the optimal basis 

SAVE separation of this vector from the preceding subbasis 

IOUT index of the critical vector (also, a fortiori, index of a discarded 
one) 

SMIN separation of the critical vector 

SEPAR separation of a discarded vector in a dependent set 

G(K,K) error in the dependent parameter a K corresponding to a unit error 
in the discarded parameter of the dependent set 

NSETC index of a dependent set after union has been performed 

DET value of the determinant of the normalized Gram matrix after removal 

of discarded vectors 

TIME value of the computation time 

WR array used to write the names of parameters involved in different 

comments issued by the subroutine and is obtained as output of the 
auxiliary subroutine TRADUC 

Auxiliary Variables 

G contains the elements of the normalized Gram matrix during the 

successive transformations into an upper triangle; G is identical to 
C at the beginning 

CLEAR logical array: CLEAR (i) = .TRUE, if the vector has already been 

used in the optimal basis (or discarded) 

IBASIS IBASIS(i) is the index of the i th vector of the basis 
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KG dimension of a subbasis 

SET logical array representing the different dependent subset: 

SET (i, j ) = -TRUE, indicates that the parameter ai belongs to the 
dependent set number j 

IR logical array: IR(i) = .TRUE, indicates that a^ is irrelevant 

Other variables are used as running subscripts or transitory storage and 
there is no need to catalog them here. 


Auxiliary Subroutine TRADUC (FOUND, Nl, K) 

This subroutine is used to translate logical information into parameter 
names in order to ease the reading of the comments. It is a part of the 
formatting of the output. The names of the parameters are given through the 
COMMON/WRITE/, in the array WA. The output of the subroutine is found in the 
array WR of the same COMMON. 


General Organization of the Subroutine COR 

COR has been divided for convenience into seven sections or logical 
units . 

Section 1- This section computes the norms AD, the normalized Gram 
matrix C and the vector (D T ) (ER) that is originally stored in the array DDA. 

It also chooses the first basis vector as being the closest to ER (the corre- 
sponding parameter would thus give the best fit if all the others had to be 
discarded). This choice will be ignored, however, if MNS(N1 + 1) is not zero. 
Control will be returned to the calling program at this point if MODE = -1. 

Section 2 - After some formatting the matrix C can be written optionally 
and initialization of the variables performed. 

Section 3- The reduction of the G matrix to an upper triangle is 
performed along with the transformation of the vector DDA. Transformation of 
columns of the matrix C is optional. At each step, the diagonal elements of 
G are searched for the biggest and the smallest. The latter is compared to 
the square of the threshold (S = THR^) . In case of dependence, control is 
transferred to section 4. Otherwise, basic and critical parameters are 
printed out with their separations. The pivotal row and the corresponding 
component of DDA are normalized (optionally the corresponding row of C) and 
the next step is started. When all the vectors have been used (IREMN = 0) 
control is transferred to section 5. 

Section 4 - The dependent parameter is discarded (MN = 0) and control 
returns to section 3 if MODE = 0. Otherwise a backward elimination is per- 
formed to obtain the components of the discarded vector on the subbasis and 
the dependent subset stored in SET. Control is returned to section 3. 
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Section 5- The solution vector DDA is finally obtained by backward 
elimination (optionally the inverse of C) . 

Section 6- When dependent sets have been found * they are printed out with 
their separations and the error in dependent parameters is computed and 
printed. Then the dependent sets are united eventually (if M0DE.GT.3). 

Section 7- In this last section some logical and formatting manipulations 
are performed and final comments are delivered. 

Remark : The analysis of the G matrix and the inversions are all 

performed in a double precision arithmetic. 

A listing of the program is given in the following pages. 
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SUBROUTINE COR (Nl,N3jTHR t.MODf 1J 
COMMON/ BAS IS/ I BAS IS ( 36 ) , P I V ( 36 ) ,G ( 36 , 36 ) 

C OMM ON/NVRT /AD( 3 6 ) , DD A ( 36 ) , C ( 36 , 36 ) 

COMMON/WRITE/WR (36) ,WA(36) 

C OMM ON/PRM/MN( 37 ),MN5(37),SENSIV(36),A(36),AMIN(36),AMAX<36) 
1,DA (36) t ER( 2400) ,DER (2400, 36) 

DIMENSION SE PAR ( 18),SET< 18? 18) , CLEAR ( 36) , CL EARB ( 36) , I R ( 36 ) 
DOUBLE PRECISION AD , DD A , G , C , GR AD, S , SA VE , SM AX , S M IN ,D E T 
LOGICAL *1 CLEAR, CLEARR, SET,IR, I NV 
DATA PEO,C 0MMA/3H 1 =,1H,/ 

CALL CLOCK( ISTART) 

MODE = MODE! 

C ^SECTION 1* COMPUTATION OF THE GRAM MATRIX 

SMAX = O.DO 
DO 500 1=1 ,N1 
AD ( I ) = 0 • DO 

IF (MN( I ) • EQ • 0 ) GO TO 2500 

1503 DO 503 K = 1 , N 3 

503 A D ( I ) = AD ( I ) + DBLE ( DER ( K, I ) )**2 

IF ( A D ( I).LT.l.D-70) MN ( I ) =0 

A D ( I ) = DSORT ( AD ( I ) ) 

2500 DO 500 J = 1 , N 1 

500 C ( I , J ) = 0 • DO 

1501 DO 501 1 = 1 T N 1 

I F ( MN ( I ) • EO * 0 ) GO TO 50 1 

1502 DO 502 J= 1 , I 

I F ( I .EQ. J) GO TO 502 

I F ( MN ( J) . EQ .0 ) GO TO 502 

1504 DO 504 K = 1 ,N3 

504 C(ItJ) = C(I,J) + DBLE ( DER ( K, I ) )*DBLE( DER ( K, J ) ) 

C ( I , J ) = C (I , J ) / ( AD ( I )*AD( J ) ) 

C( J, I ) = C ( I , J ) 

C DEFINE THE FIRST BASIC PARAMETER 

I F ( SMAX, GT. DABS (C ( I , J) ) ) GO TO 502 

SMAX = D AB S ( C ( I , J) ) 

NEXT = I 
502 CONT INUE 

C ( I , I ) = 1 • DO 
GRAD = 0 • DO 
1510 DO 510 K= 1 , N3 

510 GRAD = GRAD + DB L E ( DER ( K , I ) ) *DB L E ( ER ( K ) ) 

DDA ( I ) = GR AD /AD ( I ) 

501 CONTINUE 

N2 = N 1 + 1 

IF(MNS(N2 ) .NE.O) NEXT = MNS(N2) 


PRM 
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C ^SECTION 2* INITIALISATION 

IF(MODE.FO.-l) RETURN 

INV = MODEl.GF.lO 

IF ( INV) M DOE = M DDF 1 — 10 

MN(N2) = 1 

IF(MODE.EO.O) GO TO 507 

I F ( N1 • GT • 1 5 ) GO TO 507 

WRITE( 6,6505 ) 

F ORMA T { 1H1 , 10X, ' PARAMETERS CORRELATION MATRIX*) 

00 506 1=1 ,N1 

WR I TE ( 6 , 6506 ) WA(I),(C(I,J),J=1,I) 

FORMAT( /1X,A4,2X, 15F8.4) 

IF(MODE.LT.A) GO TO 1508 

W R I T E ( 6 , 60 1 ) W A ( NEXT) 

F ORMAT ( 1 1 BASIC »,30X,* CRITICAL 1 /* PARAMET ER * , 5 X , • S EPAR AT I ON • 

1 , 15X, * PARAMETER » ,5X, • SE PAR AT I ON* /4X, A4/) 

C INITIALISE 

1508 DO 508 1=1, Ml 

1509 DO 509 J = 1 , N 1 

IF(I .LE.18.AND. J.LE.18) SET(I,J) = .FALSE. 

G(I,J) = C ( I , J ) 

509 C ( I , J) = 0 • DO 
C< 1,1) = 1 • DO 
I R ( I ) = MN ( I ) . EO . 0 
MN ( I ) = ( CMN( I ) + l ) /2 )*2 

CL EAR ( I ) = .FALSE. 

IF ( .NOT. IR ( I ) ) GO TO 508 

DDA ( I ) = 0 • DO 
CL b AR ( I ) = .TRUE. 

508 MN ( N2 ) = M I MO ( M N ( N 2 ) , MN { I ) ) 

DET = 1 • DO 
S = THR *#2 
KG = 0 
NSET = 0 
P I V < NE XT ) = 1. 


505 

6505 
1506 

506 

6506 

507 

601 
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C * SECT I ON 3* ANALYSIS OF. THE VECTOR SET 

C THE NEW BASIS VECTOR IS NEXT 

1 CLEAR (NEXT ) = .TRUE. 

IF ( I R ( NEXT ) ) GOTO 2004 

KG = KG + 1 
IBASIS(KG) = NEXT 
IF (KG.GE .NI ) GO TO 20 

C COMPUTE THE NEW ARRAY OF REMNANT VECTORS 

1002 DO 2 J= 1 » N 1 

I F ( CLEAR ( J ) ) GO TO 2 

C TRANSFORM INPUT VECTOR DDA 

DDA(J) = DDA(J) - DDA( NEXT ) *G ( J, NEX T ) 

I F ( .NOT. INV) GO TO 1003 

1203 DO 203 1=1, Nl 

203 C(J,I) = C(J,I) - C(NEXT, I )*G( J,NEXT) 

1003 DO 3 1=1, J 

IF (CLEAR ( I ) ) GO TO 3 

G(J,I) = G ( J , I ) - G(J,NEXT)*G(NEXT» I ) 

G ( I , J ) = G ( J , I ) 

3 CONTINUE 

2 CONTINUE 

C FIND THE DEPENDENT, THE OPTIMAL AND THE CRITICAL VECTORS 
2004 IREMN = 0 

SMAX = O.DO 
SMIN = l.DO 

1004 DO 4 1=1, Nl 

IF (CLEAR ( I ) ) GO TO 4 

IREMN = IRFMN + 1 
SAVE = G ( I , I ) 

IF (SAVE. GT. SMIN) GO TO 5 

C STORE THE CRITICAL 
SMIN = SAVE 
I OUT = I 

5 IF ( SAVE. LE. SMAX ) GO TO 4 

C STORE THE UPTIMAL 
SMAX = SAVE 
NEXT = I 

4 CONTINUE 

I F ( IREMN. EO.O) GO TO 20 

C CHECK THE DEPENDENT 

IF(SMIN.LE.S) GOTO 8 

C COMPUTE THE VALUE OF THE DETERMINANT 
PIV(NEXT) = SMAX 
DET = DET*SMA>C 

C NORMALISE ROW AND COLUMN NEXT BEFORE THE NEW CYCLE 
SAVE = DSORT ( SMAX ) 

IF (M0DE.LT.4 ) GO TO 1007 

SMIN = DSORT (SMIN) 

WRI TE ( 6,600 ) WA ( NEXT ) , SAVE,WA( I OUT ) , SMIN 
600 F0RMATI4X , A4, 8 X , 1 PE 9 . 2 , 1 8 X , A4 , 8X.E9.2) 

1007 DO 7 1=1, Nl 

IF(INV) C ( NEXT , I ) = CtNEXT, I ) /SMAX 

IF (CLEAR ( I ) ) GO TO 7 

G ( NEXT , I ) = G(NEXT ,1 )/ SMAX 
7 CONTINUE 

DDA (NEXT ) = DDA ( NEXT ) /SMAX 

GO TO 1 
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C ^SECTION 4* ANALYSIS OF THE DEPENDENCE WITHIN THE SUB-BASIS 

8 MSET = NSET + 1 

DDA ( I OUT ) = O.DO 

MN( I OUT ) = 0 

CL EAR ( I OUT ) = .TRUE. 

IF (MODE.EQ.O) GO TO 2004 

SET ( I OUT, NSET) = .TRUE. 

C FIND THE COMPONENTS OF TOUT ON THE BASIS 

1009 DO 9 1 = 1 1 KG 
KBACK = KG+l-I 

K = I BAS I S < KBACK ) 

IF< I .EO. 1 ) GO TO 9 

1010 DO 10 J= 2 , I 
KBACK1 = KG- I + J 

K1 = I BAS I S ( KBACK1 ) 

10 G ( K , I OUT ) = G ( K » I OU T ) - G ( K, K 1 ) * G < K 1 , I OUT ) 

9 CONTINUE 

1011 DO 11 1=1, KG 
K = I B A S I S ( I ) 

GRAD = SMIN + (G(K, I0UT)**2)*PI V ( K ) 

I F ( DABS(GRAD) ♦ L E • S ) GO TO 11 
C MEMORISE THE DEPENDENT VECTOR IN THE LOGICAL ARRAY ''SET' 1 
SET ( K , NSET ) = .TRUE. 

C SET THE MAGIC NUMBER TO 1 FOR DEPENDENCE OF THE NEXT BASIS VECTOR WITH I0U1 
MN< K ) = 1 

11 CONTINUE 

C MEMORISE THE SEPARATION OF THE DEPENDENT VECTOR IN THE ARRAY •■SEPAR 1 ' 


SEPAR(NSET) = D SGR T ( DMAX 1 ( 0 . DO, S MI N > ) 

GO TO 2004 

C ^SECTION 5* SOLVE THE EQUATION G*DDA = D'*ER < INVERT G EVENTIJALLY)- 

20 IF(KG.EO.l) GO TO 1121 


1120 DO 120 1=2, KG 
KBACK = KG+l-I 

K = IBASIS(KBACK) 

2120 DO 120 J = 2, I 

KBACK1 = KG- I + J 
K 1 = I BA S I S ( KBACK 1 ) 

IF ( .NOT. INV) GO TO 120 

1220 DO 220 LI = 1,KG 
L = IBASIS(Ll) 

220 C ( K , L ) = C ( K , L ) - G ( K , K 1 ) *C ( K 1 , L ) 

120 DDA(K) = DDA(K) - G ( K , K 1 ) *DDA ( K 1 ) 

1121 DO 121 1 = 1, KG 
K = I BAS IS( I ) 

121 DDA(K) = DD A ( K ) / AD ( K ) 
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C * SECT I ON 6* FIND THE STRUCTURE OF THE SET 

IF(MODE.EO.O) GO TO 9999 

IF(NSET.EQ.O) GO TO 24 

IFIMODE.LT. 5) GO TO 24 

C WRITE THE DEPENDENT SETS 
WRITE! 6,621 ) 

621 F0RMAT(/11X, 'DEPENDENT SETS OF PARAMETERS’/ 1 X , • SE PARA T I ON ' / ) 

1022 DO 22 J= 1 » NSET 

1023 DO 23 1=1, N1 

IF(MN( I ) .EO.O.AND.SET! I , J ) ) IOUT = I 

23 CLEAR! I ) = SET! I , J ) 

CALL TRADUC (CLEAR ,N1,NW) 

WRITE! 6,622) SEPAR(J), ( WR ( I ) , C OMMA , 1= 1 ,NW ) 

622 FORMAT! 1 X , E9 . 2 / ( 1 OX , 24 ( A4 , A 1 )/)) 

IF (MODE. LT. 6) GO TO 22 

K = 0 

1012 DO 12 1=1, N1 

I F ( MN ( I) .EO.O ) GO TO 12 

IF ( .NOT. SET ( I , J ) ) GO TO 12 

K = K + 1 

G ( K, K ) = -G ( I , I OUT )*AD( I OUT) /AD! I ) 

WR (K) = W A ( I ) 

12 CONTINUE 

14 WR I TE ( 6,614) W A ( I OUT ) , PE 0 , ( WR ( I ) , PEO , G ( I , I > , I = 1 , K ) 

614 FORMAT!/' IF *,A4,A3,' 1, THEN ' , 6! 3 X, A 4, A3 , 1 PE9 . 2 ) / 

l(19X,6(3X,A4,A3,E9.2)/)> 

22 CONTINUE 

24 WRITE! 6,624) N 1 , THR, (MN ( K ) , K= 1 ,N1 ) 

624 F ORM A T ( ' 1 ***TFT£RE ARE*, 13, » PARAMETERS IN THIS PROBLEM . 

1 THE SEPARATION THRESHOLD WAS',E9.2 /• MAGIC NUMBER ',3612/) 

I F ( NSET .EO.O) GO TO 32 

IF(M0DE.LT.3) GO TO 32 

C*** ********** ********** UNION OF THE DEPENDENT SETS **************** 
1125 DO 125 1=1, Ml 
125 CLEAR ( I ) = .FALSE. 

NSETC = 0 

1025 DO 25 J= 1 « NSET 

IF ( CLEAR ( J ) ) GO TO 25 

NSETC = NSETC + 1 

1026 DO 26 K= J , NSET 

I F ( CLEAR ( K ) ) GO TO 26 

1027 DO 27 1=1, N1 

IF (SET! I , J ) . AND. SET! I ,K) ) GO TO 1028 

27 CONTINUE 

GO TO 26 

1028 DO 28 1=1, N1 

28 SET ( I , NSETC) = SET ( I , J) . OR. SET! I ,K ) 

SEP AR ( NSETC ) = AMINl! SEPAR ( J) , SEPAR (K) ) 

CLEAR ( K ) = .TRUE. 

26 CONTINUE 

25 CONTINUE 


****** 
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C ^SECTION 7* PREPARE AND WRITE THE FINAL DIAGNOSTIC 

102 9 no 29 J= 1 » NSE TC 
1030 no 30 1=1, Ml 

CLEAR ( I ) = (MN( I ).EQ. 1) .AND. SET ( I,J) 

30 CL E ARBI I ) = ( MN ( I ) . EO . 0 ) . AND. S ET ( I , J ) 

CALL TRADUC (CLEAR, N 1 , NW ) 

WRITE! 6,631 ) J,SFPAR(J),(WR(I) , COMMA, 1=1 ,NW) 

631 FORMAT!//' ***0 E PE NO ENT SET NUMBER', 13,' ***SE PA R AT I ON =',E9.2/ 

1* *•/• * A TRUE VALUE IS OBTAINED FOR '/(• * ' » 8X ,24 ( A4 , A 1 )/) ) 

CALL TRADUC (CLEARB, Ml ,NW ) 

WRITE! 6,629) (WR( I ) ,COMMA, I=1,NW) 

629 FORMAT!' * IF IS KNOWN THE TRUE VALUE OF '/(• * • , 8X , 24 ( A4 , A 1 )/) ) 

29 CONTINUE 

32 IFIMODE.LT. 2) RETURN 

CWRITE THE NAME OF INDEPENDENT, IRRELEVANT, DROPPED AND UNUSED PARAMETERS 
1033 DO 33 1=1, N1 

CLEAR! I ) = MN! I ) .E0.2 

33 CLEARB ( I ) = I R ( I ) . AND . ( M N S ( I ) . M E . 0 ) 

CALL TRADUC (CLEAR ,N1,NW ) 

WRITE(6,634) ( W R ( I ),COMMA, I=1,NW> 

634 FORMAT!//' * INDEPENDENT PARAMETERS '/(• * • , 2 X, 24! A4 , A 1 )/)) 

CALL TRADUC (CLEARB, N1 ,NW ) 

WRI TE ( 6,63 5) ( WR ( I ) , C OMMA , I = 1 , NW ) 

635 FORMAT!' *'/' * IRRELEVANT PARAMETERS '/(• * ' , 2X , 24 ( A4 , A 1 )/)) 
1036 DO 36 1=1, N1 

CL EAR ( I ) = MN(I ).EO.O.AND.MNS(I ). ME . 0 
36 CLEARB! I ) = MNS(I).E0.0 

CALL TRADUC (CLEAR, N 1 , NW ) 

WRITE! 6,637) ( WR ( I ) , C OMMA , I = 1 , NW ) 

637 FORMAT!' *'/• * NOT ESTIMATED'/!' * ' , 2X , 24 ( A4 , A 1 )/)) 

CALL TRADUC (CLEARB ,N1 ,NW ) 

WRITE(6,638) ( WR ( I ) , COMMA, 1= 1,NW) 

638 FORMAT!' *'/' * NOT USED'/!' * ' , 2X , 2 4 ( A 4, A 1 )/)) 

CALL CLOCK! I END ) 

TIME = ( IEND-ISTART)/100. 

WRITE(6,606) DET ,T IME 

606 FORMAT!///' DET= • , D22 . 15 , • TIME =',F9.3,' SEC') 

9999 RETURN 
END 
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SUBROUTINE TRADUC(FOUNO,Nl,K) 
C0MM0N/WRITE/WR(36> ,WA(36) 
DIMENSION FOUND(Nl) 

LOGICAL FOUND 

DATA BLANK, EMPTY/' ','NONE'/ 

K = 0 

1001 DO 1 1=1, N1 

WR ( I ) = BLANK 

IF ( .NOT. FOUND ( I ) ) GO TO 1 

K = K + 1 
W R ( K ) = WA( I ) 

1 CONTINUE 

IF(K.GT.O) RETURN 

K= 1 

W R ( 1 ) = EMPTY 

RETURN 

END 
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EXAMPLE OF COMPUTER OUTPUT 


To illustrate the capabilities of the program, two examples of output are 
given on figures 17 and 18. In the first example, a test matrix DER with 
known relationships between the columns was input to the subroutine COR to 
simulate a case of nonuniqueness where the model had 14 parameters (designated 
here by the symbolic names Ax, A 2 , . . ., A^) . The threshold was set equal 
to 0.01 and A 10 was discarded a priori (MNS(10) set to 0) and A 5 was found 
to be irrelevant. Figure 17 shows first the construction of the optimal basis. 
It indicates that the first parameter (chosen by the user) is A 2 - Then A 4 
is selected with a separation equal to 1 (indicating orthogonality between the 
vectors 3Y/8A 2 and atf/SAi*). At the same time A 3 is found to be critical 
with a separation of 0.347. The third parameter selected is Aq with again 
a separation of 1, and A 3 is still critical. At this point, however, A 9 
was found to be dependent and was discarded, but this is displayed later on. 
This process continues, the basis vectors becoming less and less orthogonal 
(property 0B3 of the optimal basis shown in appendix C) . 

Some parameters are missing in the first column; indeed they have been 
discarded either during the analysis (as dependent parameters) or before (as 
irrelevant or discarded a priori by the user right from the beginning) . As 
for the discarded parameters, they appear in the next display, where we are 
told that A g and A 9 are dependent with a separation of 0.0025. A 9 was dis- 
carded and its value will be unchanged. The value of Aq that will be 
obtained in the identification is a function of the value of A 9 , and the next 
line indicates that if the value of A 9 is increased by one unit (absolute 
variation - symbolized by a prime) then Aq will decrease by 40 units. This 
information is very useful to estimate 'the error in Aq because of the lack 
of knowledge of A 9 , and also it can be used to correct a posteriori the value 
of A 8 if a better value has been obtained for A 9 from other measurements, 
without running the whole identification again. 

Then another dependent set is found with three parameters. A*, A 2 , and 
A 3 . The parameter A 3 was discarded and the next line indicates how a varia- 
tion in A 3 will affect the values found for Ai and A 2 . The last dependent 
set includes four parameters and Am was discarded. Note that Aq appears 
again, therefore the first and the last set will be united in the final 
results . 

The final results are given next on this output. Two dependent sets 
appear; A 9 and Ax 4 have been discarded from the first and their values should 
be known in order to obtain correct values for Aq, Ah, and Ax 3 . In the same 
way A 3 was discarded from the second set. Then the names of the independent 
and irrelevant parameters are displayed. Finally, four parameters have been 
discarded and will not be estimated in this identification, and also, the 
computer reminds us that Axo was not used at all because of the a priori 
decision to do so at the beginning of the identification. The last piece of 
information concerns the value of the determinant of the reduced Gram matrix 
(that is after suppression of rows and columns corresponding to the discarded 
parameters) and, gives the computation time in seconds. All of these results 
agreed exactly with the known initial data. 
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Another example is shown on figure 18. It corresponds to the real case 
(case A) of the platform identification, when the reaction wheels are not used 
(natural damping). The dependence between the parameters Jxz, Jx, Jz and a ^3 
is found and the dependence coefficients are denoted by the "prime" following 
the name, so that the last line of the first set of comments should read: if 

Jxz is incremented by one unit, then the value obtained for Jx will change 
by an amount of 2.55 10~ l unit 3 Jz by 4.0 units and a\ 3 will be decreased 
by 4.85 10~ & unit. Note that the analysis of this 13 parameter problem only 
required 2 sec of computation time. 
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TABLE 1.- RESULTS OF THE IDENTIFICATION FOR TWO DIFFERENT 


IDENTIFICATION RUNS IN CASE A (FREE OSCILLATION CASE) 


Name 

Value 

Error bounds 

Units 

Sensitivity 

Dependency 

index 


Run 

1: Final results when no 

(threshold 

parameter is discarded 
= 0) 

Jx 

1. 

23 

10 2 

±1 . 6 



slug 

-ft 2 

12.4 

2 

J y 

1 

22 

10 2 

±2.7 





6.7 

2 

Jz 

1 

60 

10 2 

±19 





1.3 

2 

Jxy 

-3 

24 


±1.4 





.32 

2 

Jxz 

9. 

64 


±5 





.54 

2 

Jyz 

DMP X 

-1. 

3. 

31 

25 

10‘ 2 

±5.4 

±2.2 

10' 

2 

ft- 11 

3-sec 

.036 

.037 

2 

2 

DMPy 

3. 

09 

IO -2 

±2.8 

10" 

2 



.018 

2 

DMP, 

1. 

19 

IO' 1 

±6.9 

io- 

2 



.016 

2 

by 

9. 

87 

10“ 3 

±1.8 

io- 

3 

ft 

-lb 

.057 

2 

bz 

7. 

77 

10” 1 

±7.6 

io- 

3 



12.5 

2 

a I 3 

-2. 

53 

io- 4 

±2.8 

io- 

5 



.13 

2 

Run . 

2: 

Final results when 

^xz i- s 

discarded by the 

dependence 





analysis 

(threshold = 0.02) 


Jx 

1. 

18 

10 2 

±1.3 



slug 

-ft 2 

12.7 

1 

Jy 

1 . 

20 

10 2 

±1.8 





7.7 

2 

Jz 

1 . 

21 

10 2 

±7.9 





1.06 

1 

Jxy 

- 1 . 

70 


±1.2 





. 17 

2 

^xz 

- 1 . 

67 


0 





.091 

0 

^YZ 

2. 

44 

io- 2 

±2.6 





.00085 

2 

DMP X 

2. 

78 

10~ 2 

±1 . 6 

io- 

2 

ft- lb-sec 

.032 

2 

DMPy 

5. 

14 

IO" 2 

±2.5 

io- 

2 

i 


.030 

2 

DMPv 

8. 

03 

IO -2 

±3.5 

10" 

2 

i 


.014 

2 

b y 

1 . 

00 

IO" 2 

±1.5 

io- 

3 

ft 

-lb 

.055 

2 

bz 

7. 

73 

10“ 1 

±6.8 

10" 

3 



12.5 

2 

a l 3 

- 1 . 

89 

IO" 4 

±1.4 

10" 

5 



.15 

1 
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TABLE 2.- RESULTS OF THE IDENTIFICATION IN CASE B (REACTION 

WHEELS DAMPING) 


Name 

Value 

Error bounds 

Units 

Sensitivity 

Jx 

119 

±5.3 

slug-ft 2 

1.32 

J y 

125 

±5 



1.7 

Jz 

213 

±130 



.02 

^xy 

4.7 

±2.4 



.2 

T 

J xz 

-8.0 

±2.7 



.03 

^yz 

-6.1 

±2.6 



.02 

“i 

67 

±4 

ft- 11 

3-sec 

1.4 

a 2 

121 

±3.5 



.7 

“3 

147 

±16 



.13 

b y 

-3.4 10“ 3 

±1. 10" 3 

ft 

-lb 

.09 

b z 

-.76 

±1.3 10” 2 



3.7 

*13 

-1.0 10" 3 

±3. 10" 4 



.2 

a l 6 

.137 

±2. 10" 2 



.38 

a l 7 

. 135 

1+ 

I—* 

o 

1 

CjO 



1.2 

a l 8 

. 122 

0 



.0 

Cl 

4.24 10" 2 

±1.3 10" 2 

cps 

.13 

C2 

3.30 10" 2 

±1.3 10- 2 



. 13 

c 3 

4.40 10- 2 

±1.8 10" 2 

’ 

\ 

.14 

a 22 

1.36 

±2. 10' 2 

ft- 11 

D-sec 

4.5 

a 2 3 

1.35 

±3. 10“ 2 



4.0 

a 24 

1.00 

0 



0 


Dependency 

index 

2 


2 

2 

2 

2 

0 
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Figure 7 .- 


Linear closeness in 


two tensions. 
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G = D T D 

Grom determinant 

DET(G) = V 2 


Recursive computations of 


V 


2 


V k 2 + , = V 2 x S 2 + , 
DET k+1 = DET k x Si. 



DET(G) = DET k x 


Sii = S 


2 

k + I 




First basis vectors 


a 


Figure 11.- Computation technique. 




Figure 12.- General arrangement for the identification of the satellite simulator. 
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Roll rate, deg/sec 


.80 


Run 2, S 0 =0.02 



Time, sec 


Figure 15.- Measured and computed time histories of the roll rate of the platform at the end of two 

different identification runs (Case A) . 




Time, sec 


Figure 16.- Measured and computed time histories of the roll rate (Case B) 






BASIC 



CRITICAL 


PARAMETER 

SEPARATION 


PARAMETER 

SEPARATION 

A?. 





A4 

l.OOE 00 


A3 

3.47E-01 

A8 

l.OOE 00 


A3 

3.47E-01 

A12 

l.OOE 00 


A3 

3.47E-01 

Al 

9.57E-01 


A3 

3.47E-01 

A6 

7 .11 E-01 


A 14 

4.66E-01 

All 

7.07E-0 1 


A 7 

9.90E-02 

A 1 3 

7.07E-01 


A7 

9.90E-02 

A7 

9.90E-02 


A 7 

9.90E-02 


DEPENDENT SETS OF PARAMETERS 



SFPARAT I ON 





0.25E-02 

Aft, A 9 , 




IF A9 » = 

1, THEN A8 ’ =-4.01E 

01 



0.38E-02 

Al, A2 , A3, 




ii 

m 

< 

LL 

1, THEN Al' =-2.01 E 

00 

A 2 1 = -l.OOE 

01 

0.71E-08 

Aft, All, A 1 3, Al 4, 




IF A 1 4 • = 

1, THFN Aft* =-2 . OOE 

00 

All' =-3 .OOF 

00 A 1 3* 


444THERE ARE 14 PARAMETERS TN THIS PROBLEM . 
MAGIC NIJMBFR 11020 221001210 

THE SEPARATION THRFSHOLO WAS 0.10F-01 


444DE PENDENT SET NUMBER 1 **4SE PAR AT I ON = 0.71F-0B 

4 

* A TRUE VALUE IS OBTAINED FOR 

* Aft, All, A IB , 

* IF IS KNOWN THF TRUE VALUE OF 

4 A9 , A 14 , 


444DF PENDENT SET NUMBER 2 *44SE PAR AT I ON = 0.3BF-0? 

4 A TRUE VALUE IS OBTAINED FOR 
4 Al, A2 , 

4 IF IS KNOWN THF TRUE VALUE OF 
4 A3, 


4 INDEPENDENT PARAMETERS 
4 A4, A6, A 1, A12, 

4 IRRELEVANT PARAMETERS 
4 A5 , 

4 

* NOT ESTIMATED 
4 A3, A5 , A9, A 14, 

4 

4 NOT USED 
4 A10, 

D FT = 0. 1 13448R44BR4488D-02 TIME = 0.734 SEC 


Figure 17.- Computer analysis in the test case of 14 parameters. 
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PARAMETER 
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PARAMETER 

SEPARATION 
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BY 

l.OOD 00 

JXZ 

2.77D-01 

DMPY 

9.77D-01 

JXZ 

2. 510-01 

JY 

9.23D-01 

JXZ 

2.510-01 

DMPZ 

8.07D-01 

JYZ 

1.02D-01 

DM PX 

5 .380-01 

JYZ 

9 .570-02 

JX 

5.16D-01 

JYZ 

9.55D-02 

A 13 

4.55D-01 

JXZ 

3 .040-02 

JYZ 

4.27D-02 

JXY 

3.130-02 

BZ 

3. 11D-02 

JXY 

2.79D-02 

JXY 

2.77D-02 

JXY 

2.77D-02 


DEPENDENT SETS OF PARAMETERS 



SEPARATION 




0.48E-02 

JX, JZ, JXZ, A 1 3 , 



IF JXZ' = 

1, THEN JX* = 2.55D-01 

JZ 1 = 4.00 D 

00 A 1 3 » 


***THERE ARE 
MAGIC NUMBER 


13 PARAMETERS IN THIS PROBLEM 
1212022220221 


SEPARATION THRESHOLD WAS 0.20E- 


01 


***DEPENDENT SET NUMBER 1 ***SE PAR ATI ON = 0.48E-0? 

* A TRUE VALUE IS OBTAINED FOR 

* JX, JZ, A13, 

* IF IS KNOWN THE TRUE VALUE OF 

* JXZ , 


* INDEPENDENT PARAMETERS 

* JY, JXY, JYZ,DMPX,DMPY,DMPZ, BY, BZ, 

* 
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* N ONE , 

* NOT ESTIMATED 

* JXZ, 

* NOT USED 

* BX, 


DFT= 0.113607241258916D-10 TIME 


1.910 SEC 


Figure 18.- Computer analysis of the identification of the parameters of the 
platform in the free oscillation case (Case A, run 2 ). 


74 


NASA-Langley, 1971 19 A- 3 80 7 


II II 


ii ilium i mi in 


mi ■ ii 


ii i ii 


i 


ii 


i 



